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

Merge branch 'develop'

This commit is contained in:
Atsushi Togo 2025-06-12 16:22:01 +09:00
parent f0c26cec8a dd60053a74
commit bbd055d9b7
52 changed files with 3141 additions and 2216 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.pre-commit-config.yaml
View File

@ -11,7 +11,7 @@ repos:
exclude: ^example/AlN-LDA/
- repo: https://github.com/astral-sh/ruff-pre-commit
rev: v0.11.7
rev: v0.11.13
hooks:
- id: ruff
args: [ "--fix", "--show-fixes" ]

75
c/_phono3py.cpp
View File

@ -50,18 +50,17 @@ static Darray *convert_to_darray(nb::ndarray<> npyary) {
// printf("Data shift:%lu [%lu, %lu]\n", adrs_shift, i_sigma, i_temp);
// }
void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
nb::ndarray<> py_g_zero, nb::ndarray<> py_frequencies,
nb::ndarray<> py_eigenvectors,
nb::ndarray<> py_triplets,
nb::ndarray<> py_bz_grid_addresses,
nb::ndarray<> py_D_diag, nb::ndarray<> py_Q,
nb::ndarray<> py_fc3, nb::ndarray<> py_svecs,
nb::ndarray<> py_multi, nb::ndarray<> py_masses,
nb::ndarray<> py_p2s_map, nb::ndarray<> py_s2p_map,
nb::ndarray<> py_band_indices, int64_t symmetrize_fc3_q,
int64_t make_r0_average, nb::ndarray<> py_all_shortest,
double cutoff_frequency, int64_t openmp_per_triplets) {
void py_get_interaction(
nb::ndarray<> py_fc3_normal_squared, nb::ndarray<> py_g_zero,
nb::ndarray<> py_frequencies, nb::ndarray<> py_eigenvectors,
nb::ndarray<> py_triplets, nb::ndarray<> py_bz_grid_addresses,
nb::ndarray<> py_D_diag, nb::ndarray<> py_Q, nb::ndarray<> py_fc3,
nb::ndarray<> py_fc3_nonzero_indices, nb::ndarray<> py_svecs,
nb::ndarray<> py_multi, nb::ndarray<> py_masses, nb::ndarray<> py_p2s_map,
nb::ndarray<> py_s2p_map, nb::ndarray<> py_band_indices,
int64_t symmetrize_fc3_q, int64_t make_r0_average,
nb::ndarray<> py_all_shortest, double cutoff_frequency,
int64_t openmp_per_triplets) {
Darray *fc3_normal_squared;
Darray *freqs;
_lapack_complex_double *eigvecs;
@ -72,6 +71,7 @@ void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
int64_t *D_diag;
int64_t (*Q)[3];
double *fc3;
char *fc3_nonzero_indices;
double (*svecs)[3];
int64_t (*multi)[2];
double *masses;
@ -100,6 +100,7 @@ void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
} else {
is_compact_fc3 = 1;
}
fc3_nonzero_indices = (char *)py_fc3_nonzero_indices.data();
svecs = (double (*)[3])py_svecs.data();
for (i = 0; i < 2; i++) {
multi_dims[i] = py_multi.shape(i);
@ -113,9 +114,10 @@ void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
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,
all_shortest, cutoff_frequency, openmp_per_triplets);
fc3_nonzero_indices, is_compact_fc3, svecs,
multi_dims, multi, masses, p2s, s2p, band_indices,
symmetrize_fc3_q, make_r0_average, all_shortest,
cutoff_frequency, openmp_per_triplets);
free(fc3_normal_squared);
fc3_normal_squared = NULL;
@ -129,8 +131,9 @@ void py_get_pp_collision(
nb::ndarray<> py_triplets, nb::ndarray<> py_triplet_weights,
nb::ndarray<> py_bz_grid_addresses, nb::ndarray<> py_bz_map,
int64_t bz_grid_type, nb::ndarray<> py_D_diag, nb::ndarray<> py_Q,
nb::ndarray<> py_fc3, nb::ndarray<> py_svecs, nb::ndarray<> py_multi,
nb::ndarray<> py_masses, nb::ndarray<> py_p2s_map, nb::ndarray<> py_s2p_map,
nb::ndarray<> py_fc3, nb::ndarray<> py_fc3_nonzero_indices,
nb::ndarray<> py_svecs, nb::ndarray<> py_multi, nb::ndarray<> py_masses,
nb::ndarray<> py_p2s_map, nb::ndarray<> py_s2p_map,
nb::ndarray<> py_band_indices, nb::ndarray<> py_temperatures_THz,
int64_t is_NU, int64_t symmetrize_fc3_q, int64_t make_r0_average,
nb::ndarray<> py_all_shortest, double cutoff_frequency,
@ -147,6 +150,7 @@ void py_get_pp_collision(
int64_t *D_diag;
int64_t (*Q)[3];
double *fc3;
char *fc3_nonzero_indices;
double (*svecs)[3];
int64_t (*multi)[2];
double *masses;
@ -176,6 +180,7 @@ void py_get_pp_collision(
} else {
is_compact_fc3 = 1;
}
fc3_nonzero_indices = (char *)py_fc3_nonzero_indices.data();
svecs = (double (*)[3])py_svecs.data();
for (i = 0; i < 2; i++) {
multi_dims[i] = py_multi.shape(i);
@ -191,9 +196,10 @@ void py_get_pp_collision(
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_THz, is_NU, symmetrize_fc3_q,
make_r0_average, all_shortest, cutoff_frequency, openmp_per_triplets);
D_diag, Q, fc3, fc3_nonzero_indices, is_compact_fc3, svecs, multi_dims,
multi, masses, p2s, s2p, band_indices, temperatures_THz, is_NU,
symmetrize_fc3_q, make_r0_average, all_shortest, cutoff_frequency,
openmp_per_triplets);
free(band_indices);
band_indices = NULL;
@ -206,7 +212,8 @@ void py_get_pp_collision_with_sigma(
nb::ndarray<> py_frequencies, nb::ndarray<> py_eigenvectors,
nb::ndarray<> py_triplets, nb::ndarray<> py_triplet_weights,
nb::ndarray<> py_bz_grid_addresses, nb::ndarray<> py_D_diag,
nb::ndarray<> py_Q, nb::ndarray<> py_fc3, nb::ndarray<> py_svecs,
nb::ndarray<> py_Q, nb::ndarray<> py_fc3,
nb::ndarray<> py_fc3_nonzero_indices, nb::ndarray<> py_svecs,
nb::ndarray<> py_multi, nb::ndarray<> py_masses, nb::ndarray<> py_p2s_map,
nb::ndarray<> py_s2p_map, nb::ndarray<> py_band_indices,
nb::ndarray<> py_temperatures_THz, int64_t is_NU, int64_t symmetrize_fc3_q,
@ -222,6 +229,7 @@ void py_get_pp_collision_with_sigma(
int64_t *D_diag;
int64_t (*Q)[3];
double *fc3;
char *fc3_nonzero_indices;
double (*svecs)[3];
int64_t (*multi)[2];
double *masses;
@ -249,6 +257,7 @@ void py_get_pp_collision_with_sigma(
} else {
is_compact_fc3 = 1;
}
fc3_nonzero_indices = (char *)py_fc3_nonzero_indices.data();
svecs = (double (*)[3])py_svecs.data();
for (i = 0; i < 2; i++) {
multi_dims[i] = py_multi.shape(i);
@ -264,8 +273,8 @@ void py_get_pp_collision_with_sigma(
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_THz, is_NU, symmetrize_fc3_q,
fc3_nonzero_indices, is_compact_fc3, svecs, multi_dims, multi, masses,
p2s, s2p, band_indices, temperatures_THz, is_NU, symmetrize_fc3_q,
make_r0_average, all_shortest, cutoff_frequency, openmp_per_triplets);
free(band_indices);
@ -795,13 +804,11 @@ int64_t py_tpl_get_triplets_reciprocal_mesh_at_q(
return num_ir;
}
int64_t py_tpl_get_BZ_triplets_at_q(nb::ndarray<> py_triplets,
int64_t grid_point,
nb::ndarray<> py_bz_grid_address,
nb::ndarray<> py_bz_map,
nb::ndarray<> py_map_triplets,
nb::ndarray<> py_D_diag, nb::ndarray<> py_Q,
int64_t bz_grid_type) {
int64_t py_tpl_get_BZ_triplets_at_q(
nb::ndarray<> py_triplets, int64_t grid_point,
nb::ndarray<> py_bz_grid_address, nb::ndarray<> py_bz_map,
nb::ndarray<> py_map_triplets, nb::ndarray<> py_D_diag, nb::ndarray<> py_Q,
nb::ndarray<> py_reciprocal_lattice, int64_t bz_grid_type) {
int64_t (*triplets)[3];
int64_t (*bz_grid_address)[3];
int64_t *bz_map;
@ -809,6 +816,7 @@ int64_t py_tpl_get_BZ_triplets_at_q(nb::ndarray<> py_triplets,
int64_t num_map_triplets;
int64_t *D_diag;
int64_t (*Q)[3];
double (*reciprocal_lattice)[3];
int64_t num_ir;
triplets = (int64_t (*)[3])py_triplets.data();
@ -818,10 +826,11 @@ int64_t py_tpl_get_BZ_triplets_at_q(nb::ndarray<> py_triplets,
num_map_triplets = (int64_t)py_map_triplets.shape(0);
D_diag = (int64_t *)py_D_diag.data();
Q = (int64_t (*)[3])py_Q.data();
reciprocal_lattice = (double (*)[3])py_reciprocal_lattice.data();
num_ir = ph3py_get_BZ_triplets_at_q(triplets, grid_point, bz_grid_address,
bz_map, map_triplets, num_map_triplets,
D_diag, Q, bz_grid_type);
num_ir = ph3py_get_BZ_triplets_at_q(
triplets, grid_point, bz_grid_address, bz_map, map_triplets,
num_map_triplets, D_diag, Q, reciprocal_lattice, bz_grid_type);
return num_ir;
}

7
c/gridsys.c
View File

@ -138,7 +138,7 @@ int64_t gridsys_get_bz_grid_addresses(
if (!niggli_reduce(niggli_lattice, GRIDSYS_NIGGLI_TOLERANCE)) {
return 0;
}
if (!lagmat_inverse_matrix_d3(inv_Lr, (double(*)[3])niggli_lattice,
if (!lagmat_inverse_matrix_d3(inv_Lr, (double (*)[3])niggli_lattice,
GRIDSYS_NIGGLI_TOLERANCE)) {
return 0;
}
@ -156,7 +156,7 @@ int64_t gridsys_get_bz_grid_addresses(
bzgrid->bzg2grg = bzg2grg;
bzgrid->type = bz_grid_type;
lagmat_multiply_matrix_l3(bzgrid->Q, inv_Mpr_int, Q);
lagmat_copy_matrix_d3(bzgrid->reclat, (double(*)[3])niggli_lattice);
lagmat_copy_matrix_d3(bzgrid->reclat, (double (*)[3])niggli_lattice);
for (i = 0; i < 3; i++) {
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = PS[i];
@ -220,7 +220,7 @@ int64_t gridsys_get_bz_triplets_at_q(
const int64_t (*bz_grid_addresses)[3], const int64_t *bz_map,
const int64_t *map_triplets, const int64_t num_map_triplets,
const int64_t D_diag[3], const int64_t Q[3][3],
const int64_t bz_grid_type) {
const double reciprocal_lattice[3][3], const int64_t bz_grid_type) {
RecgridConstBZGrid *bzgrid;
int64_t i, j, num_ir;
@ -237,6 +237,7 @@ int64_t gridsys_get_bz_triplets_at_q(
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
for (j = 0; j < 3; j++) {
bzgrid->reclat[i][j] = reciprocal_lattice[i][j];
bzgrid->Q[i][j] = Q[i][j];
}
}

3
c/gridsys.h
View File

@ -272,7 +272,8 @@ int64_t gridsys_get_bz_triplets_at_q(
int64_t (*ir_triplets)[3], const int64_t bz_grid_index,
const int64_t (*bz_grid_addresses)[3], const int64_t *bz_map,
const int64_t *map_triplets, const int64_t num_map_triplets,
const int64_t D_diag[3], const int64_t Q[3][3], const int64_t bz_grid_type);
const int64_t D_diag[3], const int64_t Q[3][3],
const double reciprocal_lattice[3][3], const int64_t bz_grid_type);
/**
* @brief Return integration weight of linear tetrahedron method

37
c/phono3py.c
View File

@ -63,12 +63,13 @@ int64_t ph3py_get_interaction(
const _lapack_complex_double *eigenvectors, const int64_t (*triplets)[3],
const int64_t num_triplets, const int64_t (*bz_grid_addresses)[3],
const int64_t D_diag[3], const int64_t Q[3][3], const double *fc3,
const int64_t is_compact_fc3, const double (*svecs)[3],
const int64_t multi_dims[2], const int64_t (*multiplicity)[2],
const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
const int64_t *band_indices, const int64_t symmetrize_fc3_q,
const int64_t make_r0_average, const char *all_shortest,
const double cutoff_frequency, const int64_t openmp_per_triplets) {
const char *fc3_nonzero_indices, const int64_t is_compact_fc3,
const double (*svecs)[3], const int64_t multi_dims[2],
const int64_t (*multiplicity)[2], const double *masses,
const int64_t *p2s_map, const int64_t *s2p_map, const int64_t *band_indices,
const int64_t symmetrize_fc3_q, const int64_t make_r0_average,
const char *all_shortest, const double cutoff_frequency,
const int64_t openmp_per_triplets) {
RecgridConstBZGrid *bzgrid;
AtomTriplets *atom_triplets;
int64_t i, j;
@ -102,6 +103,7 @@ int64_t ph3py_get_interaction(
atom_triplets->s2p_map = s2p_map;
atom_triplets->make_r0_average = make_r0_average;
atom_triplets->all_shortest = all_shortest;
atom_triplets->nonzero_indices = fc3_nonzero_indices;
itr_get_interaction(fc3_normal_squared, g_zero, frequencies,
(lapack_complex_double *)eigenvectors, triplets,
@ -127,7 +129,8 @@ int64_t ph3py_get_pp_collision(
const int64_t (*bz_grid_addresses)[3], /* thm */
const int64_t *bz_map, /* thm */
const int64_t bz_grid_type, const int64_t D_diag[3], const int64_t Q[3][3],
const double *fc3, const int64_t is_compact_fc3, const double (*svecs)[3],
const double *fc3, const char *fc3_nonzero_indices,
const int64_t is_compact_fc3, const double (*svecs)[3],
const int64_t multi_dims[2], const int64_t (*multiplicity)[2],
const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
const Larray *band_indices, const Darray *temperatures_THz,
@ -169,6 +172,7 @@ int64_t ph3py_get_pp_collision(
atom_triplets->s2p_map = s2p_map;
atom_triplets->make_r0_average = make_r0_average;
atom_triplets->all_shortest = all_shortest;
atom_triplets->nonzero_indices = fc3_nonzero_indices;
ppc_get_pp_collision(imag_self_energy, relative_grid_address, frequencies,
(lapack_complex_double *)eigenvectors, triplets,
@ -192,13 +196,14 @@ int64_t ph3py_get_pp_collision_with_sigma(
const int64_t (*triplets)[3], const int64_t num_triplets,
const int64_t *triplet_weights, const int64_t (*bz_grid_addresses)[3],
const int64_t D_diag[3], const int64_t Q[3][3], const double *fc3,
const int64_t is_compact_fc3, const double (*svecs)[3],
const int64_t multi_dims[2], const int64_t (*multiplicity)[2],
const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
const Larray *band_indices, const Darray *temperatures_THz,
const int64_t is_NU, const int64_t symmetrize_fc3_q,
const int64_t make_r0_average, const char *all_shortest,
const double cutoff_frequency, const int64_t openmp_per_triplets) {
const char *fc3_nonzero_indices, const int64_t is_compact_fc3,
const double (*svecs)[3], const int64_t multi_dims[2],
const int64_t (*multiplicity)[2], const double *masses,
const int64_t *p2s_map, const int64_t *s2p_map, const Larray *band_indices,
const Darray *temperatures_THz, const int64_t is_NU,
const int64_t symmetrize_fc3_q, const int64_t make_r0_average,
const char *all_shortest, const double cutoff_frequency,
const int64_t openmp_per_triplets) {
RecgridConstBZGrid *bzgrid;
AtomTriplets *atom_triplets;
int64_t i, j;
@ -232,6 +237,7 @@ int64_t ph3py_get_pp_collision_with_sigma(
atom_triplets->s2p_map = s2p_map;
atom_triplets->make_r0_average = make_r0_average;
atom_triplets->all_shortest = all_shortest;
atom_triplets->nonzero_indices = fc3_nonzero_indices;
ppc_get_pp_collision_with_sigma(
imag_self_energy, sigma, sigma_cutoff, frequencies,
@ -405,7 +411,7 @@ int64_t ph3py_get_BZ_triplets_at_q(
const int64_t (*bz_grid_addresses)[3], const int64_t *bz_map,
const int64_t *map_triplets, const int64_t num_map_triplets,
const int64_t D_diag[3], const int64_t Q[3][3],
const int64_t bz_grid_type) {
const double reciprocal_lattice[3][3], const int64_t bz_grid_type) {
RecgridConstBZGrid *bzgrid;
int64_t i, j, num_ir;
@ -423,6 +429,7 @@ int64_t ph3py_get_BZ_triplets_at_q(
bzgrid->PS[i] = 0;
for (j = 0; j < 3; j++) {
bzgrid->Q[i][j] = Q[i][j];
bzgrid->reclat[i][j] = reciprocal_lattice[i][j];
}
}
bzgrid->size = num_map_triplets;

34
c/phono3py.h
View File

@ -53,12 +53,13 @@ int64_t ph3py_get_interaction(
const _lapack_complex_double *eigenvectors, const int64_t (*triplets)[3],
const int64_t num_triplets, const int64_t (*bz_grid_addresses)[3],
const int64_t D_diag[3], const int64_t Q[3][3], const double *fc3,
const int64_t is_compact_fc3, const double (*svecs)[3],
const int64_t multi_dims[2], const int64_t (*multi)[2],
const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
const int64_t *band_indices, const int64_t symmetrize_fc3_q,
const int64_t make_r0_average, const char *all_shortest,
const double cutoff_frequency, const int64_t openmp_per_triplets);
const char *fc3_nonzero_indices, const int64_t is_compact_fc3,
const double (*svecs)[3], const int64_t multi_dims[2],
const int64_t (*multi)[2], const double *masses, const int64_t *p2s_map,
const int64_t *s2p_map, const int64_t *band_indices,
const int64_t symmetrize_fc3_q, const int64_t make_r0_average,
const char *all_shortest, const double cutoff_frequency,
const int64_t openmp_per_triplets);
int64_t ph3py_get_pp_collision(
double *imag_self_energy,
const int64_t relative_grid_address[24][4][3], /* thm */
@ -68,7 +69,8 @@ int64_t ph3py_get_pp_collision(
const int64_t (*bz_grid_addresses)[3], /* thm */
const int64_t *bz_map, /* thm */
const int64_t bz_grid_type, const int64_t D_diag[3], const int64_t Q[3][3],
const double *fc3, const int64_t is_compact_fc3, const double (*svecs)[3],
const double *fc3, const char *fc3_nonzero_indices,
const int64_t is_compact_fc3, const double (*svecs)[3],
const int64_t multi_dims[2], const int64_t (*multi)[2],
const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
const Larray *band_indices, const Darray *temperatures_THz,
@ -81,13 +83,14 @@ int64_t ph3py_get_pp_collision_with_sigma(
const int64_t (*triplets)[3], const int64_t num_triplets,
const int64_t *triplet_weights, const int64_t (*bz_grid_addresses)[3],
const int64_t D_diag[3], const int64_t Q[3][3], const double *fc3,
const int64_t is_compact_fc3, const double (*svecs)[3],
const int64_t multi_dims[2], const int64_t (*multi)[2],
const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
const Larray *band_indices, const Darray *temperatures_THz,
const int64_t is_NU, const int64_t symmetrize_fc3_q,
const int64_t make_r0_average, const char *all_shortest,
const double cutoff_frequency, const int64_t openmp_per_triplets);
const char *fc3_nonzero_indices, const int64_t is_compact_fc3,
const double (*svecs)[3], const int64_t multi_dims[2],
const int64_t (*multi)[2], const double *masses, const int64_t *p2s_map,
const int64_t *s2p_map, const Larray *band_indices,
const Darray *temperatures_THz, const int64_t is_NU,
const int64_t symmetrize_fc3_q, const int64_t make_r0_average,
const char *all_shortest, const double cutoff_frequency,
const int64_t 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 int64_t (*triplets)[3],
@ -170,7 +173,8 @@ int64_t ph3py_get_BZ_triplets_at_q(
int64_t (*triplets)[3], const int64_t grid_point,
const int64_t (*bz_grid_addresses)[3], const int64_t *bz_map,
const int64_t *map_triplets, const int64_t num_map_triplets,
const int64_t D_diag[3], const int64_t Q[3][3], const int64_t bz_grid_type);
const int64_t D_diag[3], const int64_t Q[3][3],
const double reciprocal_lattice[3][3], const int64_t bz_grid_type);
int64_t ph3py_get_integration_weight(
double *iw, char *iw_zero, const double *frequency_points,
const int64_t num_band0, const int64_t relative_grid_address[24][4][3],

19
c/real_to_reciprocal.c
View File

@ -75,7 +75,7 @@ static lapack_complex_double get_pre_phase_factor(
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] */
/* fc3_reciprocal[num_patom, 3, num_patom, 3, num_patom, 3] */
void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
const double q_vecs[3][3], const double *fc3,
const int64_t is_compact_fc3,
@ -277,7 +277,7 @@ static void real_to_reciprocal_elements(
const int64_t pi0, const int64_t pi1, const int64_t pi2,
const int64_t leg_index) {
int64_t i, j, k, l;
int64_t num_satom, adrs_shift;
int64_t num_satom, adrs_shift_atoms, adrs_shift_fc3;
lapack_complex_double phase_factor;
double fc3_rec_real[27], fc3_rec_imag[27];
@ -309,8 +309,11 @@ static void real_to_reciprocal_elements(
continue;
}
}
adrs_shift =
i * 27 * num_satom * num_satom + j * 27 * num_satom + k * 27;
adrs_shift_atoms = i * num_satom * num_satom + j * num_satom + k;
if (!atom_triplets->nonzero_indices[adrs_shift_atoms]) {
continue;
}
adrs_shift_fc3 = adrs_shift_atoms * 27;
phase_factor =
phonoc_complex_prod(phase_factor1[j], phase_factor2[k]);
@ -320,19 +323,19 @@ static void real_to_reciprocal_elements(
for (l = 0; l < 27; l++) {
fc3_rec_real[l] +=
lapack_complex_double_real(phase_factor) *
fc3[adrs_shift + l] * 3;
fc3[adrs_shift_fc3 + l] * 3;
fc3_rec_imag[l] +=
lapack_complex_double_imag(phase_factor) *
fc3[adrs_shift + l] * 3;
fc3[adrs_shift_fc3 + l] * 3;
}
} else {
for (l = 0; l < 27; l++) {
fc3_rec_real[l] +=
lapack_complex_double_real(phase_factor) *
fc3[adrs_shift + l];
fc3[adrs_shift_fc3 + l];
fc3_rec_imag[l] +=
lapack_complex_double_imag(phase_factor) *
fc3[adrs_shift + l];
fc3[adrs_shift_fc3 + l];
}
}
}

1
c/real_to_reciprocal.h
View File

@ -48,6 +48,7 @@ typedef struct {
const int64_t *s2p_map;
int64_t make_r0_average;
const char *all_shortest;
const char *nonzero_indices; // for compact fc3
} AtomTriplets;
void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,

15
c/recgrid.c
View File

@ -193,20 +193,11 @@ double recgrid_get_tolerance_for_BZ_reduction(const RecgridBZGrid *bzgrid) {
int64_t i, j;
double tolerance;
double length[3];
double reclatQ[3][3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
reclatQ[i][j] = bzgrid->reclat[i][0] * bzgrid->Q[0][j] +
bzgrid->reclat[i][1] * bzgrid->Q[1][j] +
bzgrid->reclat[i][2] * bzgrid->Q[2][j];
}
}
for (i = 0; i < 3; i++) {
length[i] = 0;
for (j = 0; j < 3; j++) {
length[i] += reclatQ[j][i] * reclatQ[j][i];
length[i] += bzgrid->reclat[j][i] * bzgrid->reclat[j][i];
}
length[i] /= bzgrid->D_diag[i] * bzgrid->D_diag[i];
}
@ -233,8 +224,8 @@ RecgridMats *recgrid_alloc_RotMats(const int64_t size) {
rotmats->size = size;
if (size > 0) {
if ((rotmats->mat = (int64_t(*)[3][3])malloc(sizeof(int64_t[3][3]) *
size)) == NULL) {
if ((rotmats->mat = (int64_t (*)[3][3])malloc(sizeof(int64_t[3][3]) *
size)) == NULL) {
warning_print("Memory could not be allocated ");
warning_print("(RecgridMats, line %d, %s).\n", __LINE__, __FILE__);
free(rotmats);

328
c/reciprocal_to_normal.c
View File

@ -49,16 +49,31 @@
#include "lapack_wrapper.h"
#ifdef MEASURE_R2N
#include <time.h>
#include <unistd.h>
#endif
static void reciprocal_to_normal_squared_no_threading(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *e0, const lapack_complex_double *e1,
const lapack_complex_double *e2, const int64_t *band_indices,
const int64_t num_band, const double cutoff_frequency);
static void reciprocal_to_normal_squared_g_threading(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *e0, const lapack_complex_double *e1,
const lapack_complex_double *e2, const int64_t *band_indices,
const int64_t num_band, const double cutoff_frequency);
static double get_fc3_sum(const lapack_complex_double *e0,
const lapack_complex_double *e1,
static double get_fc3_sum(const lapack_complex_double *e1,
const lapack_complex_double *e2,
const lapack_complex_double *fc3_reciprocal,
const int64_t num_band);
static double get_fc3_sum_blas_like(const lapack_complex_double *e0,
const lapack_complex_double *e1,
const lapack_complex_double *e2,
const lapack_complex_double *fc3_reciprocal,
const int64_t num_band);
// Testing efficiency of BLAS
#ifdef MULTITHREADED_BLAS
static double get_fc3_sum_blas(const lapack_complex_double *e0,
const lapack_complex_double *e1,
@ -66,11 +81,7 @@ static double get_fc3_sum_blas(const lapack_complex_double *e0,
const lapack_complex_double *fc3_reciprocal,
const int64_t num_band);
#endif
static double get_fc3_sum_blas_like(const lapack_complex_double *e0,
const lapack_complex_double *e1,
const lapack_complex_double *e2,
const lapack_complex_double *fc3_reciprocal,
const int64_t num_band);
void reciprocal_to_normal_squared(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
@ -80,15 +91,16 @@ void reciprocal_to_normal_squared(
const lapack_complex_double *eigvecs2, const double *masses,
const int64_t *band_indices, const int64_t num_band,
const double cutoff_frequency, const int64_t openmp_per_triplets) {
int64_t i, j, ij, num_atom, use_multithreaded_blas;
int64_t i, j, ij, num_atom;
double *inv_sqrt_masses;
lapack_complex_double *e0, *e1, *e2;
/* Inverse sqrt mass is multiplied with eigenvectors to reduce number
* of */
/* operations in get_fc3_sum. Three eigenvector matrices are looped
* by */
/* first loop leveraging contiguous memory layout of [e0, e1, e2].
/* Inverse sqrt mass is multiplied with eigenvectors to reduce
* number of */
/* operations in get_fc3_sum. Three eigenvector matrices are
* looped by */
/* first loop leveraging contiguous memory layout of [e0,
* e1, e2].
*/
num_atom = num_band / 3;
inv_sqrt_masses = (double *)malloc(sizeof(double) * num_band);
@ -132,64 +144,17 @@ void reciprocal_to_normal_squared(
free(inv_sqrt_masses);
inv_sqrt_masses = NULL;
#ifdef MEASURE_R2N
double loopTotalCPUTime, loopTotalWallTime;
time_t loopStartWallTime;
clock_t loopStartCPUTime;
#endif
#ifdef MEASURE_R2N
loopStartWallTime = time(NULL);
loopStartCPUTime = clock();
#endif
#ifdef MULTITHREADED_BLAS
if (openmp_per_triplets) {
use_multithreaded_blas = 0;
reciprocal_to_normal_squared_no_threading(
fc3_normal_squared, g_pos, num_g_pos, fc3_reciprocal, freqs0,
freqs1, freqs2, e0, e1, e2, band_indices, num_band,
cutoff_frequency);
} else {
use_multithreaded_blas = 1;
reciprocal_to_normal_squared_g_threading(
fc3_normal_squared, g_pos, num_g_pos, fc3_reciprocal, freqs0,
freqs1, freqs2, e0, e1, e2, band_indices, num_band,
cutoff_frequency);
}
#else
use_multithreaded_blas = 0;
#ifdef _OPENMP
#pragma omp parallel for if (!openmp_per_triplets)
#endif
#endif
for (i = 0; i < num_g_pos; i++) {
if (freqs0[band_indices[g_pos[i][0]]] > cutoff_frequency &&
freqs1[g_pos[i][1]] > cutoff_frequency &&
freqs2[g_pos[i][2]] > cutoff_frequency) {
#ifdef MULTITHREADED_BLAS
if (use_multithreaded_blas) {
fc3_normal_squared[g_pos[i][3]] =
get_fc3_sum_blas(e0 + band_indices[g_pos[i][0]] * num_band,
e1 + g_pos[i][1] * num_band,
e2 + g_pos[i][2] * num_band,
fc3_reciprocal, num_band) /
(freqs0[band_indices[g_pos[i][0]]] * freqs1[g_pos[i][1]] *
freqs2[g_pos[i][2]]);
} else {
#endif
fc3_normal_squared[g_pos[i][3]] =
get_fc3_sum_blas_like(
e0 + band_indices[g_pos[i][0]] * num_band,
e1 + g_pos[i][1] * num_band,
e2 + g_pos[i][2] * num_band, fc3_reciprocal, num_band) /
(freqs0[band_indices[g_pos[i][0]]] * freqs1[g_pos[i][1]] *
freqs2[g_pos[i][2]]);
#ifdef MULTITHREADED_BLAS
}
#endif
} else {
fc3_normal_squared[g_pos[i][3]] = 0;
}
}
#ifdef MEASURE_R2N
loopTotalCPUTime = (double)(clock() - loopStartCPUTime) / CLOCKS_PER_SEC;
loopTotalWallTime = difftime(time(NULL), loopStartWallTime);
printf(" %1.3fs (%1.3fs CPU)\n", loopTotalWallTime, loopTotalCPUTime);
#endif
free(e0);
e0 = NULL;
@ -197,41 +162,177 @@ void reciprocal_to_normal_squared(
e2 = NULL;
}
static double get_fc3_sum(const lapack_complex_double *e0,
const lapack_complex_double *e1,
// This function is more efficient than the one with multithreading
// getting rid of no concurrency over g.
static void reciprocal_to_normal_squared_no_threading(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *e0, const lapack_complex_double *e1,
const lapack_complex_double *e2, const int64_t *band_indices,
const int64_t num_band, const double cutoff_frequency) {
int64_t i, j, k, ll, bi_prev;
lapack_complex_double *fc3_e0, fc3_elem, zero;
zero = lapack_make_complex_double(0, 0);
bi_prev = -1;
fc3_e0 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) *
num_band * num_band);
for (i = 0; i < num_g_pos; i++) {
if (freqs0[band_indices[g_pos[i][0]]] > cutoff_frequency &&
freqs1[g_pos[i][1]] > cutoff_frequency &&
freqs2[g_pos[i][2]] > cutoff_frequency) {
if (bi_prev != g_pos[i][0]) {
bi_prev = g_pos[i][0];
for (j = 0; j < num_band * num_band; j++) {
fc3_e0[j] = zero;
}
for (j = 0; j < num_band; j++) {
for (k = 0; k < num_band; k++) {
for (ll = 0; ll < num_band; ll++) {
fc3_elem = phonoc_complex_prod(
fc3_reciprocal[j * num_band * num_band +
k * num_band + ll],
e0[band_indices[g_pos[i][0]] * num_band + j]);
fc3_e0[k * num_band + ll] =
lapack_make_complex_double(
lapack_complex_double_real(
fc3_e0[k * num_band + ll]) +
lapack_complex_double_real(fc3_elem),
lapack_complex_double_imag(
fc3_e0[k * num_band + ll]) +
lapack_complex_double_imag(fc3_elem));
}
}
}
}
fc3_normal_squared[g_pos[i][3]] =
get_fc3_sum(e1 + g_pos[i][1] * num_band,
e2 + g_pos[i][2] * num_band, fc3_e0, num_band) /
(freqs0[band_indices[g_pos[i][0]]] * freqs1[g_pos[i][1]] *
freqs2[g_pos[i][2]]);
} else {
fc3_normal_squared[g_pos[i][3]] = 0;
}
}
free(fc3_e0);
fc3_e0 = NULL;
}
// This is less efficient than the one without multithreading
// but can be called when unit cell is large.
static void reciprocal_to_normal_squared_g_threading(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *e0, const lapack_complex_double *e1,
const lapack_complex_double *e2, const int64_t *band_indices,
const int64_t num_band, const double cutoff_frequency) {
int64_t i;
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_g_pos; i++) {
if (freqs0[band_indices[g_pos[i][0]]] > cutoff_frequency &&
freqs1[g_pos[i][1]] > cutoff_frequency &&
freqs2[g_pos[i][2]] > cutoff_frequency) {
fc3_normal_squared[g_pos[i][3]] =
get_fc3_sum_blas_like(e0 + band_indices[g_pos[i][0]] * num_band,
e1 + g_pos[i][1] * num_band,
e2 + g_pos[i][2] * num_band,
fc3_reciprocal, num_band) /
(freqs0[band_indices[g_pos[i][0]]] * freqs1[g_pos[i][1]] *
freqs2[g_pos[i][2]]);
} else {
fc3_normal_squared[g_pos[i][3]] = 0;
}
}
}
static double get_fc3_sum(const lapack_complex_double *e1,
const lapack_complex_double *e2,
const lapack_complex_double *fc3_reciprocal,
const lapack_complex_double *fc3_e0,
const int64_t num_band) {
int64_t i, j, jk;
int64_t i, j;
double sum_real, sum_imag;
lapack_complex_double e_012_fc3, fc3_i_e_12, *e_12_cache;
const lapack_complex_double *fc3_i;
lapack_complex_double *fc3_e0_e1, fc3_elem;
e_12_cache = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) *
num_band * num_band);
fc3_e0_e1 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) *
num_band);
for (i = 0; i < num_band; i++) {
fc3_e0_e1[i] = lapack_make_complex_double(0, 0);
}
for (i = 0; i < num_band; i++) {
for (j = 0; j < num_band; j++) {
fc3_elem = phonoc_complex_prod(fc3_e0[i * num_band + j], e1[i]);
fc3_e0_e1[j] = lapack_make_complex_double(
lapack_complex_double_real(fc3_e0_e1[j]) +
lapack_complex_double_real(fc3_elem),
lapack_complex_double_imag(fc3_e0_e1[j]) +
lapack_complex_double_imag(fc3_elem));
}
}
sum_real = 0;
sum_imag = 0;
for (i = 0; i < num_band; i++) {
for (j = 0; j < num_band; j++) {
e_12_cache[i * num_band + j] = phonoc_complex_prod(e1[i], e2[j]);
}
fc3_elem = phonoc_complex_prod(fc3_e0_e1[i], e2[i]);
sum_real += lapack_complex_double_real(fc3_elem);
sum_imag += lapack_complex_double_imag(fc3_elem);
}
free(fc3_e0_e1);
fc3_e0_e1 = NULL;
return (sum_real * sum_real + sum_imag * sum_imag);
}
static double get_fc3_sum_blas_like(const lapack_complex_double *e0,
const lapack_complex_double *e1,
const lapack_complex_double *e2,
const lapack_complex_double *fc3_reciprocal,
const int64_t num_band) {
int64_t i, j;
double sum_real, sum_imag, retval_real, retval_imag;
lapack_complex_double *e_12, fc3_e_12, fc3_e_012;
e_12 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) *
num_band * num_band);
for (i = 0; i < num_band; i++) {
fc3_i = fc3_reciprocal + i * num_band * num_band;
for (jk = 0; jk < num_band * num_band; jk++) {
fc3_i_e_12 = phonoc_complex_prod(fc3_i[jk], e_12_cache[jk]);
e_012_fc3 = phonoc_complex_prod(e0[i], fc3_i_e_12);
sum_real += lapack_complex_double_real(e_012_fc3);
sum_imag += lapack_complex_double_imag(e_012_fc3);
memcpy(e_12 + i * num_band, e2,
sizeof(lapack_complex_double) * num_band);
for (j = 0; j < num_band; j++) {
e_12[i * num_band + j] =
phonoc_complex_prod(e1[i], e_12[i * num_band + j]);
}
}
free(e_12_cache);
e_12_cache = NULL;
return (sum_real * sum_real + sum_imag * sum_imag);
retval_real = 0;
retval_imag = 0;
for (i = 0; i < num_band; i++) {
sum_real = 0;
sum_imag = 0;
for (j = 0; j < num_band * num_band; j++) {
fc3_e_12 = phonoc_complex_prod(
fc3_reciprocal[i * num_band * num_band + j], e_12[j]);
sum_real += lapack_complex_double_real(fc3_e_12);
sum_imag += lapack_complex_double_imag(fc3_e_12);
}
fc3_e_012 = phonoc_complex_prod(
e0[i], lapack_make_complex_double(sum_real, sum_imag));
retval_real += lapack_complex_double_real(fc3_e_012);
retval_imag += lapack_complex_double_imag(fc3_e_012);
}
free(e_12);
e_12 = NULL;
return retval_real * retval_real + retval_imag * retval_imag;
}
#ifdef MULTITHREADED_BLAS
@ -271,46 +372,3 @@ static double get_fc3_sum_blas(const lapack_complex_double *e0,
lapack_complex_double_imag(retval);
}
#endif
static double get_fc3_sum_blas_like(const lapack_complex_double *e0,
const lapack_complex_double *e1,
const lapack_complex_double *e2,
const lapack_complex_double *fc3_reciprocal,
const int64_t num_band) {
int64_t i, j;
double sum_real, sum_imag, retval_real, retval_imag;
lapack_complex_double *e_12, fc3_e_12, fc3_e_012;
e_12 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) *
num_band * num_band);
for (i = 0; i < num_band; i++) {
memcpy(e_12 + i * num_band, e2, 16 * num_band);
for (j = 0; j < num_band; j++) {
e_12[i * num_band + j] =
phonoc_complex_prod(e1[i], e_12[i * num_band + j]);
}
}
retval_real = 0;
retval_imag = 0;
for (i = 0; i < num_band; i++) {
sum_real = 0;
sum_imag = 0;
for (j = 0; j < num_band * num_band; j++) {
fc3_e_12 = phonoc_complex_prod(
fc3_reciprocal[i * num_band * num_band + j], e_12[j]);
sum_real += lapack_complex_double_real(fc3_e_12);
sum_imag += lapack_complex_double_imag(fc3_e_12);
}
fc3_e_012 = phonoc_complex_prod(
e0[i], lapack_make_complex_double(sum_real, sum_imag));
retval_real += lapack_complex_double_real(fc3_e_012);
retval_imag += lapack_complex_double_imag(fc3_e_012);
}
free(e_12);
e_12 = NULL;
return retval_real * retval_real + retval_imag * retval_imag;
}

12
c/triplet_grid.c
View File

@ -262,7 +262,7 @@ static void get_BZ_triplets_at_q_type1(int64_t (*triplets)[3],
int64_t bzgp[3], G[3];
int64_t bz_adrs0[3], bz_adrs1[3], bz_adrs2[3];
const int64_t *gp_map;
const int64_t(*bz_adrs)[3];
const int64_t (*bz_adrs)[3];
double d2, min_d2, tolerance;
double LQD_inv[3][3];
@ -348,13 +348,14 @@ static void get_BZ_triplets_at_q_type2(int64_t (*triplets)[3],
int64_t bzgp[3], G[3];
int64_t bz_adrs0[3], bz_adrs1[3], bz_adrs2[3];
const int64_t *gp_map;
const int64_t(*bz_adrs)[3];
const int64_t (*bz_adrs)[3];
double d2, min_d2, tolerance;
double LQD_inv[3][3];
gp_map = bzgrid->gp_map;
bz_adrs = bzgrid->addresses;
get_LQD_inv(LQD_inv, bzgrid);
/* This tolerance is used to be consistent to BZ reduction in bzgrid. */
tolerance = recgrid_get_tolerance_for_BZ_reduction((RecgridBZGrid *)bzgrid);
@ -423,10 +424,15 @@ static void get_LQD_inv(double LQD_inv[3][3],
int64_t i, j, k;
/* LQD^-1 */
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
LQD_inv[i][j] = 0;
}
}
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
LQD_inv[i][k] =
LQD_inv[i][k] +=
bzgrid->reclat[i][j] * bzgrid->Q[j][k] / bzgrid->D_diag[k];
}
}

7
doc/changelog.md
View File

@ -2,6 +2,10 @@
# Change Log
## Jun-12-2025: Version 3.16.0
- Release to follow the change of phonopy
## Apr-30-2025: Version 3.15.1
- Release to follow the change of phonopy
@ -357,8 +361,7 @@ loader.
`FORCES_FC3` is assumed.
- TURBOMOLE interface is provided by Antti Karttunen (`--turbomole`).
- Compatibility of `fc2.hdf5` and `force_constants.hdf5` was improved for all
calculators to store physical unit information in the hdf5 file. See
{ref}`file_format_compatibility`.
calculators to store physical unit information in the hdf5 file.
## Mar-24-2019: Version 1.16.0

144
doc/command-options.md
View File

@ -197,7 +197,6 @@ created from `FORCES_FC2` and `phono3py_disp.yaml` instead of `FORCES_FC3` and
```
(sp_option)=
### `--sp` or `--save-params`
Instead of `FORCES_FC3`, `phono3py_params.yaml` is generated. This option must
@ -329,7 +328,6 @@ information about primitive cell (`primitive_matrix` key) in
```
(random_displacements_option)=
### `--rd` (`RANDOM_DISPLACEMENTS`), `--rd-fc2` (`RANDOM_DISPLACEMENTS_FC2`) and `--random-seed` (`RANDOM_SEED`)
**`phono3py-load` doesn't have this option.**
@ -596,7 +594,7 @@ expected.
### `--sigma` (`SIGMA`)
$\sigma$ value of Gaussian function for smearing when calculating imaginary part
of self energy. See the detail at {ref}`brillouinzone_sum`.
of self energy.
Multiple $\sigma$ values are also specified by space separated numerical values.
This is used when we want to test several $\sigma$ values simultaneously.
@ -806,48 +804,8 @@ $\Gamma_\lambda(\omega_\lambda) =
\Gamma^\text{U}_\lambda(\omega_\lambda)$
and this is used to calculate thermal conductivity in single-mode RTA. The
separation, i.e., the choice of G-vector, is made based on the first Brillouin
zone.
zone. See {ref}`iofile_kappa_hdf5_gamma_NU`.
The data are stored in `kappa-mxxx(-gx-sx-sdx).hdf5` file and accessed by
`gamma_N` and `gamma_U` keys. The shape of the arrays is the same as that of
`gamma` (see {ref}`kappa_hdf5_file_gamma`). An example (Si-PBEsol) is shown
below:
```bash
% phono3py-load --mesh 11 11 11 --fc3 --fc2 --br --nu
...
% ipython
In [1]: import h5py
In [2]: f = h5py.File("kappa-m111111.hdf5", 'r')
In [3]: list(f)
Out[3]:
['frequency',
'gamma',
'gamma_N',
'gamma_U',
'group_velocity',
'gv_by_gv',
'heat_capacity',
'kappa',
'kappa_unit_conversion',
'mesh',
'mode_kappa',
'qpoint',
'temperature',
'weight']
In [4]: f['gamma'].shape
Out[4]: (101, 56, 6)
In [5]: f['gamma_N'].shape
Out[5]: (101, 56, 6)
In [6]: f['gamma_U'].shape
Out[6]: (101, 56, 6)
```
### `--scattering-event-class` (`SCATTERING_EVENT_CLASS`)
@ -1202,68 +1160,9 @@ not found, `kappa-mxxx-gx-bx(-sx-sdx).hdf5` files for band indices are searched.
### `--write-gamma-detail` (`WRITE_GAMMA_DETAIL = .TRUE.`)
Each q-point triplet contribution to imaginary part of self energy is written
into `gamma_detail-mxxx-gx(-sx-sdx).hdf5` file. Be careful that this is large
data.
into `gamma_detail-mxxx-gx(-sx-sdx).hdf5` file. Be careful that this can be a
large file. See {ref}`iofile_gamma_detail_hdf5`.
In the output file in hdf5, following keys are used to extract the detailed
information.
```{table}
| dataset | Array shape |
| --------------------------- | -------------------------------------------------------------------------------------------------------------------------------------- |
| gamma_detail for `--ise` | (temperature, sampling frequency point, symmetry reduced set of triplets at a grid point, band1, band2, band3) in THz (without $2\pi$) |
| gamma_detail for `--br` | (temperature, symmetry reduced set of triplets at a grid point, band1, band2, band3) in THz (without $2\pi$) |
| mesh | Numbers of sampling mesh along reciprocal axes. |
| frequency_point for `--ise` | Sampling frequency points in THz (without $2\pi$), i.e., $\omega$ in $\Gamma_\lambda(\omega)$ |
| temperature | (temperature,), Temperatures in K |
| triplet | (symmetry reduced set of triplets at a grid point, 3), Triplets are given by the grid point indices (see below). |
| weight | (symmetry reduced set of triplets at a grid point,), Weight of each triplet to imaginary part of self energy |
```
Imaginary part of self energy (linewidth/2) is recovered by the following
script:
```python
import h5py
import numpy as np
gd = h5py.File("gamma_detail-mxxx-gx.hdf5")
temp_index = 30 # index of temperature
temperature = gd['temperature'][temp_index]
gamma_tp = gd['gamma_detail'][:].sum(axis=-1).sum(axis=-1)
weight = gd['weight'][:]
gamma = np.dot(weight, gamma_tp[temp_index])
```
For example, for `--br`, this `gamma` gives $\Gamma_\lambda(\omega_\lambda)$ of
the band indices at the grid point indicated by $\lambda$ at the temperature of
index 30. If any bands are degenerated, those `gamma` in
`kappa-mxxx-gx(-sx-sdx).hdf5` or `gamma-mxxx-gx(-sx-sdx).hdf5` type file are
averaged, but the `gamma` obtained here in this way are not symmetrized. Apart
from this symmetrization, the values must be equivalent between them.
To understand each contribution of triptle to imaginary part of self energy,
reading `phonon-mxxx.hdf5` is useful (see {ref}`write_phonon_option`). For
example, phonon triplets of three phonon scatterings are obtained by
```python
import h5py
import numpy as np
gd = h5py.File("gamma_detail-mxxx-gx.hdf5", 'r')
ph = h5py.File("phonon-mxxx.hdf5", 'r')
gp1 = gd['grid_point'][()]
triplets = gd['triplet'][:] # Sets of (gp1, gp2, gp3) where gp1 is fixed
mesh = gd['mesh'][:]
grid_address = ph['grid_address'][:]
q_triplets = grid_address[triplets] / mesh.astype('double')
# Phonons of triplets[2]
phonon_tp = [(ph['frequency'][i], ph['eigenvector'][i]) for i in triplets[2]]
# Fractions of contributions of triplets at this grid point and temperature index 30
gamma_sum_over_bands = np.dot(weight, gd['gamma_detail'][30].sum(axis=-1).sum(axis=-1).sum(axis=-1))
contrib_tp = [gd['gamma_detail'][30, i].sum() / gamma_sum_over_bands for i in range(len(weight))]
np.dot(weight, contrib_tp) # is one
```
(write_phonon_option)=
@ -1271,45 +1170,12 @@ np.dot(weight, contrib_tp) # is one
Phonon frequencies, eigenvectors, and grid point addresses are stored in
`phonon-mxxx.hdf5` file. {ref}`--pa <pa_option>` and {ref}`--nac <nac_option>`
may be required depending on calculation setting.
may be required depending on calculation setting. See {ref}`iofile_phonon_hdf5`.
```bash
% phono3py-load --mesh 11 11 11 --nac --write-phoonon
```
Contents of `phonon-mxxx.hdf5` are watched by:
```bash
In [1]: import h5py
In [2]: ph = h5py.File("phonon-m111111.hdf5", 'r')
In [3]: list(ph)
Out[3]: ['eigenvector', 'frequency', 'grid_address', 'mesh']
In [4]: ph['mesh'][:]
Out[4]: array([11, 11, 11], dtype=int32)
In [5]: ph['grid_address'].shape
Out[5]: (1367, 3)
In [6]: ph['frequency'].shape
Out[6]: (1367, 6)
In [7]: ph['eigenvector'].shape
Out[7]: (1367, 6, 6)
```
The first axis of `ph['grid_address']`, `ph['frequency']`, and
`ph['eigenvector']` corresponds to the number of q-points where phonons are
calculated. Here the number of phonons may not be equal to product of mesh
numbers ($1367 \neq 11^3$). This is because all q-points on Brillouin zone
boundary are included, i.e., even if multiple q-points are translationally
equivalent, those phonons are stored separately though these phonons are
physically equivalent within the equations employed in phono3py. Here Brillouin
zone is defined by Wigner–Seitz cell of reciprocal primitive basis vectors. This
is convenient to categorize phonon triplets into Umklapp and Normal scatterings
based on the Brillouin zone.
(read_phonon_option)=

303
doc/grid.md Normal file
View File

@ -0,0 +1,303 @@
(grid)=
# Grids in reciprocal space
The regular grid can be a conventional regular grid or a generalized regular
grid. Here the conventional regular grid means that the grids are cut parallel
to the reciprocal basis vectors. In most cases, the conventional regular grid is
used. In special case, e.g., for crystals with body center tetragonal symmetry,
the generalized regular grid can be useful. In phono3py, the generalized regular
grid is defined to be cut parallel to the reciprocal basis vectors of the
conventional unit cell.
Two types of grid data structure are used in phono3py. Normal regular grid
contains translationally unique grid points (regular grid). The other grid
includes the points on Brillouin zone (BZ) boundary (BZ grid).
## `BZGrid` class instance
Grid point information in reciprocal space is stored in the `BZGrid` class. This
class instance can be easily accessed in the following way.
```python
In [1]: import phono3py
In [2]: ph3 = phono3py.load("phono3py.yaml", produce_fc=False)
In [3]: ph3.mesh_numbers = [11, 11, 11]
In [4]: ph3.grid
Out[4]: <phono3py.phonon.grid.BZGrid at 0x1070f3b60>
```
It is recommended to read docstring in `BZGrid` by
```python
In [5]: help(ph3.grid)
```
Some attributes of this class are presented below.
```python
In [6]: ph3.grid.addresses.shape
Out[6]: (1367, 3)
In [7]: ph3.grid.D_diag
Out[7]: array([11, 11, 11])
In [8]: ph3.grid.P
Out[8]:
array([[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
In [9]: ph3.grid.Q
Out[9]:
array([[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
In [10]: ph3.grid.QDinv
Out[10]:
array([[0.09090909, 0. , 0. ],
[0. , 0.09090909, 0. ],
[0. , 0. , 0.09090909]])
In [11]: ph3.grid.PS
Out[11]: array([0, 0, 0])
```
The integer array `addresses` contains grid point addresses. Every grid point
address is represented by the unique series of three integers. These addresses
are converted to the q-points in fractional coordinates as explained in the
{ref}`section below<grid_address_to_q>`.
Unless generalized regular grid is employed, the other attributes are not
important. `D_diag` is equivalent to the three integer numbers of the specified
conventional regular grid. `P` and `Q` are are the left and right unimodular
matrix after Smith normal form: $\mathrm{D}=\mathrm{PAQ}$, respectively, where
$\mathrm{A}$ is the grid matrix. `D_diag` is the three diagonal elements of the
matrix $\mathrm{D}$. The grid matrix is usually a diagonal matrix, then
$\mathrm{P}$ and $\mathrm{Q}$ are chosen as identity matrix. `QDinv` is given
by $\mathrm{Q}\mathrm{D}^{-1}$. `PS` represents half-grid-shifts (usually always
`[0, 0, 0]` in phono3py).
## Find grid point index corresponding to grid point address
Grid point index corresponding to a grid point address is obtained using the
instance method `BZGrid.get_indices_from_addresses` as follows:
```python
In [1]: import phono3py
In [2]: ph3 = phono3py.load("phono3py_disp.yaml")
In [3]: ph3.mesh_numbers = [20, 20, 20]
In [4]: ph3.grid.get_indices_from_addresses([0, 10, 10])
Out[4]: 4448
```
This index number is different between phono3py version 1.x and 2.x.
To get the number corresponding to the phono3py version 1.x,
`store_dense_gp_map=False` should be specified in `phono3py.load`,
```python
In [5]: ph3 = phono3py.load("phono3py_disp.yaml", store_dense_gp_map=False)
In [6]: ph3.mesh_numbers = [20, 20, 20]
In [7]: ph3.grid.get_indices_from_addresses([0, 10, 10])
Out[7]: 4200
```
(grid_address_to_q)=
## q-points in fractional coordinates corresponding to grid addresses
For Gamma centered regular grid, q-points in fractional coordinates
are obtained by
```python
qpoints = addresses @ QDinv.T
```
For shifted regular grid (usually unused in phono3py),
```python
qpoints = (addresses * 2 + PS) @ (QDinv.T / 2.0)
```
The grid addresses are stored in `phonon-*.hdf5`. So for conventional
Gamma-centered regular grid, those information can be used to recover the
corresponding q-points. For example,
```python
In [1]: import h5py
In [2]: f = h5py.File("phonon-m111111.hdf5")
In [3]: import numpy as np
In [8]: f['grid_address'][:] @ np.diag(1.0 / f['mesh'][:])
Out[8]:
array([[ 0. , 0. , 0. ],
[ 0.09090909, 0. , 0. ],
[ 0.18181818, 0. , 0. ],
...,
[-0.27272727, -0.09090909, -0.09090909],
[-0.18181818, -0.09090909, -0.09090909],
[-0.09090909, -0.09090909, -0.09090909]], shape=(1367, 3))
```
(grid_triplets)=
## Grid point triplets
Three grid point indices are used to represent a q-point triplet. For example
the following command generates `gamma_detail-m111111-g5.hdf5`,
```bash
% phono3py-load phono3py.yaml --gp 5 --br --mesh 11 11 11 --write-gamma-detail
```
This file contains various information:
```python
In [1]: import h5py
In [2]: f = h5py.File("gamma_detail-m111111-g5.hdf5")
In [3]: list(f)
Out[3]:
['gamma_detail',
'grid_point',
'mesh',
'temperature',
'triplet',
'triplet_all',
'version',
'weight']
In [4]: f['gamma_detail'].shape
Out[4]: (101, 146, 6, 6, 6)
```
For the detailed analysis of contributions of triplets to imaginary part of
self energy a phonon mode of the grid point, it is necessary to understand the
data structure of `triplet` and `weight`.
```python
In [5]: f['triplet'].shape
Out[5]: (146, 3)
In [6]: f['weight'].shape
Out[6]: (146,)
In [7]: f['triplet'][:10]
Out[7]:
array([[ 5, 0, 6],
[ 5, 1, 5],
[ 5, 2, 4],
[ 5, 3, 3],
[ 5, 7, 10],
[ 5, 8, 9],
[ 5, 11, 118],
[ 5, 12, 117],
[ 5, 13, 116],
[ 5, 14, 115]])
```
The second index of `gamma_detail` corresponds to the first index of `triplet`.
Three integers of each triplet are the grid point indices, which means, the grid
addresses and their q-points are recovered by
```python
In [8]: import phono3py
In [9]: ph3 = phono3py.load("phono3py.yaml", produce_fc=False)
In [10]: ph3.mesh_numbers = [11, 11, 11]
In [11]: ph3.grid.addresses[f['triplet'][:10]]
Out[11]:
array([[[ 5, 0, 0],
[ 0, 0, 0],
[-5, 0, 0]],
[[ 5, 0, 0],
[ 1, 0, 0],
[ 5, 0, 0]],
[[ 5, 0, 0],
[ 2, 0, 0],
[ 4, 0, 0]],
[[ 5, 0, 0],
[ 3, 0, 0],
[ 3, 0, 0]],
[[ 5, 0, 0],
[-4, 0, 0],
[-1, 0, 0]],
[[ 5, 0, 0],
[-3, 0, 0],
[-2, 0, 0]],
[[ 5, 0, 0],
[ 0, 1, 0],
[-5, -1, 0]],
[[ 5, 0, 0],
[ 1, 1, 0],
[ 5, -1, 0]],
[[ 5, 0, 0],
[ 2, 1, 0],
[ 4, -1, 0]],
[[ 5, 0, 0],
[ 3, 1, 0],
[ 3, -1, 0]]])
n [14]: ph3.grid.addresses[f['triplet'][:10]] @ ph3.grid.QDinv
Out[14]:
array([[[ 0.45454545, 0. , 0. ],
[ 0. , 0. , 0. ],
[-0.45454545, 0. , 0. ]],
[[ 0.45454545, 0. , 0. ],
[ 0.09090909, 0. , 0. ],
[ 0.45454545, 0. , 0. ]],
[[ 0.45454545, 0. , 0. ],
[ 0.18181818, 0. , 0. ],
[ 0.36363636, 0. , 0. ]],
[[ 0.45454545, 0. , 0. ],
[ 0.27272727, 0. , 0. ],
[ 0.27272727, 0. , 0. ]],
[[ 0.45454545, 0. , 0. ],
[-0.36363636, 0. , 0. ],
[-0.09090909, 0. , 0. ]],
[[ 0.45454545, 0. , 0. ],
[-0.27272727, 0. , 0. ],
[-0.18181818, 0. , 0. ]],
[[ 0.45454545, 0. , 0. ],
[ 0. , 0.09090909, 0. ],
[-0.45454545, -0.09090909, 0. ]],
[[ 0.45454545, 0. , 0. ],
[ 0.09090909, 0.09090909, 0. ],
[ 0.45454545, -0.09090909, 0. ]],
[[ 0.45454545, 0. , 0. ],
[ 0.18181818, 0.09090909, 0. ],
[ 0.36363636, -0.09090909, 0. ]],
[[ 0.45454545, 0. , 0. ],
[ 0.27272727, 0.09090909, 0. ],
[ 0.27272727, -0.09090909, 0. ]]])
```

280
doc/hdf5_howto.md
View File

@ -1,23 +1,18 @@
(hdf5_howto)=
# How to read the results stored in hdf5 files
# Using phono3py hdf5 files
```{contents}
:depth: 3
:local:
```
## Using `h5py` in ipython
## How to use HDF5 python library
It is assumed that `python-h5py` is installed on the computer you interactively
use. In the following, how to see the contents of `.hdf5` files in the
interactive mode of Python. The basic usage of reading `.hdf5` files using
`h5py` is found at
[here](http://docs.h5py.org/en/latest/high/dataset.html#reading-writing-data>).
In the following example, an MgO result of thermal conductivity calculation
stored in `kappa-m111111.hdf5` (see {ref}`iofile_kappa_hdf5`) is loaded and
thermal conductivity tensor at 300 K is watched.
It is assumed that `python-h5py` is installed on the computer you
interactively use. In the following, how to see the contents of
`.hdf5` files in the interactive mode of Python. The basic usage of
reading `.hdf5` files using `h5py` is found at [here](http://docs.h5py.org/en/latest/high/dataset.html#reading-writing-data>).
Usually for running interactive python, `ipython` is recommended to
use but not the plain python. In the following example, an MgO result
of thermal conductivity calculation is loaded and thermal conductivity
tensor at 300 K is watched.
```bash
```python
In [1]: import h5py
In [2]: f = h5py.File("kappa-m111111.hdf5")
@ -90,256 +85,3 @@ In [11]: g = np.where(g > 0, g, -1)
In [12]: lifetime = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0)
```
(kappa_hdf5_file)=
## Details of `kappa-*.hdf5` file
Files name, e.g. `kappa-m323220.hdf5`, is determined by some
specific options. `mxxx`, show the numbers of sampling
mesh. `sxxx` and `gxxx` appear optionally. `sxxx` gives the
smearing width in the smearing method for Brillouin zone integration
for phonon lifetime, and `gxxx` denotes the grid number. Using the
command option of `-o`, the file name can be modified slightly. For
example `-o nac` gives `kappa-m323220.nac.hdf5` to
memorize the option `--nac` was used.
Currently `kappa-*.hdf5` file (not for the specific grid points)
contains the properties shown below.
### mesh
(Versions 1.10.11 or later)
The numbers of mesh points for reciprocal space sampling along
reciprocal axes, $a^*, b^*, c^*$.
### frequency
Phonon frequencies. The physical unit is THz, where THz
is in the ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band).
(kappa_hdf5_file_gamma)=
### gamma
Imaginary part of self energy of phonon bubble diagram (phonon-phonon
scattering). The physical unit is THz, where THz is in the ordinal frequency not
the angular frequency.
The array shape for all grid-points (irreducible q-points) is
(temperature, irreducible q-point, phonon band).
The array shape for a specific grid-point is
(temperature, phonon band).
Phonon lifetime ($\tau_\lambda=1/2\Gamma_\lambda(\omega_\lambda)$) may
be estimated from `gamma`. $2\pi$ has to be multiplied with
`gamma` values in the hdf5 file to convert the unit of ordinal
frequency to angular frequency. Zeros in `gamma` values mean that
those elements were not calculated such as for three acoustic modes at
$\Gamma$ point. The below is the copy-and-paste from the
previous section to show how to obtain phonon lifetime in pico
second:
```bash
In [8]: g = f['gamma'][30]
In [9]: import numpy as np
In [10]: g = np.where(g > 0, g, -1)
In [11]: lifetime = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0)
```
### gamma_isotope
Isotope scattering of $1/2\tau^\mathrm{iso}_\lambda$.
The physical unit is same as that of gamma.
The array shape is same as that of frequency.
### group_velocity
Phonon group velocity, $\nabla_\mathbf{q}\omega_\lambda$. The
physical unit is $\text{THz}\cdot\text{Angstrom}$, where THz
is in the ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band, 3 = Cartesian coordinates).
### heat_capacity
Mode-heat-capacity defined by
$$
C_\lambda = k_\mathrm{B}
\left(\frac{\hbar\omega_\lambda}{k_\mathrm{B} T} \right)^2
\frac{\exp(\hbar\omega_\lambda/k_\mathrm{B}
T)}{[\exp(\hbar\omega_\lambda/k_\mathrm{B} T)-1]^2}.
$$
The physical unit is eV/K.
The array shape is (temperature, irreducible q-point, phonon band).
(output_kappa)=
### kappa
Thermal conductivity tensor. The physical unit is W/m-K.
The array shape is (temperature, 6 = (xx, yy, zz, yz, xz, xy)).
(output_mode_kappa)=
### mode-kappa
Thermal conductivity tensors at k-stars (${}^*\mathbf{k}$):
$$
\sum_{\mathbf{q} \in {}^*\mathbf{k}} \kappa_{\mathbf{q}j}.
$$
The sum of this over ${}^*\mathbf{k}$ corresponding to
irreducible q-points divided by number of grid points gives
$\kappa$ ({ref}`output_kappa`), e.g.,:
```python
kappa_xx_at_index_30 = mode_kappa[30, :, :, 0].sum()/ weight.sum()
```
Be careful that until version 1.12.7, mode-kappa values were divided
by number of grid points.
The physical unit is W/m-K. Each tensor element is the sum of tensor
elements on the members of ${}^*\mathbf{k}$, i.e., symmetrically
equivalent q-points by crystallographic point group and time reversal
symmetry.
The array shape is (temperature, irreducible q-point, phonon band, 6 =
(xx, yy, zz, yz, xz, xy)).
### gv_by_gv
Outer products of group velocities for k-stars
(${}^*\mathbf{k}$) for each irreducible q-point and phonon band
($j$):
$$
\sum_{\mathbf{q} \in {}^*\mathbf{k}} \mathbf{v}_{\mathbf{q}j} \otimes
\mathbf{v}_{\mathbf{q}j}.
$$
The physical unit is
$\text{THz}^2\cdot\text{Angstrom}^2$, where THz is in the
ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band, 6 = (xx, yy, zz,
yz, xz, xy)).
### q-point
Irreducible q-points in reduced coordinates.
The array shape is (irreducible q-point, 3 = reduced
coordinates in reciprocal space).
### temperature
Temperatures where thermal conductivities are calculated. The physical
unit is K.
### weight
Weights corresponding to irreducible q-points. Sum of weights equals to
the number of mesh grid points.
### ave_pp
Averaged phonon-phonon interaction in $\text{eV}^2,
$P_{\mathbf{q}j}$:
$$
P_{\mathbf{q}j} = \frac{1}{(3n_\mathrm{a})^2} \sum_{\lambda'\lambda''}
|\Phi_{\lambda\lambda'\lambda''}|^2.
$$
This is not going to be calculated in the RTA thermal coductivity
calculation mode by default. To calculate this, `--full-pp` option
has to be specified (see {ref}`full_pp_option`).
### boundary_mfp
A value specified by {ref}`boundary_mfp_option`. The physical unit is
micrometer.
When `--boundary-mfp` option is explicitly specified, its value is stored here.
### kappa_unit_conversion
This is used to convert the physical unit of lattice thermal
conductivity made of `heat_capacity`, `group_velocity`, and
`gamma`, to W/m-K. In the single mode relaxation time (SMRT) method,
a mode contribution to the lattice thermal conductivity is given by
$$
\kappa_\lambda = \frac{1}{V_0} C_\lambda \mathbf{v}_\lambda \otimes
\mathbf{v}_\lambda \tau_\lambda^{\mathrm{SMRT}}.
$$
For example, $\kappa_{\lambda,{xx}}$ is calculated by:
```bash
In [1]: import h5py
In [2]: f = h5py.File("kappa-m111111.hdf5")
In [3]: kappa_unit_conversion = f['kappa_unit_conversion'][()]
In [4]: weight = f['weight'][:]
In [5]: heat_capacity = f['heat_capacity'][:]
In [6]: gv_by_gv = f['gv_by_gv'][:]
In [7]: gamma = f['gamma'][:]
In [8]: kappa_unit_conversion * heat_capacity[30, 2, 0] * gv_by_gv[2, 0] / (2 * gamma[30, 2, 0])
Out[8]:
array([ 1.02050241e+03, 1.02050241e+03, 1.02050241e+03,
4.40486382e-15, 0.00000000e+00, -4.40486382e-15])
In [9]: f['mode_kappa'][30, 2, 0]
Out[9]:
array([ 1.02050201e+03, 1.02050201e+03, 1.02050201e+03,
4.40486209e-15, 0.00000000e+00, -4.40486209e-15])
```
## How to know grid point index number corresponding to grid address
Running with `--write-gamma`, hdf5 files are written out with file names
such as `kappa-m202020-g4448.hdf5`. You may want to know the grid point
index number with given grid address. This is done as follows:
```bash
In [1]: import phono3py
In [2]: ph3 = phono3py.load("phono3py_disp.yaml")
In [3]: ph3.mesh_numbers = [20, 20, 20]
In [4]: ph3.grid.get_indices_from_addresses([0, 10, 10])
Out[4]: 4448
```
This index number is different between phono3py version 1.x and 2.x.
To get the number corresponding to the phono3py version 1.x,
`store_dense_gp_map=False` should be specified in `phono3py.load`,
```bash
In [5]: ph3 = phono3py.load("phono3py_disp.yaml", store_dense_gp_map=False)
In [6]: ph3.mesh_numbers = [20, 20, 20]
In [7]: ph3.grid.get_indices_from_addresses([0, 10, 10])
Out[7]: 4200
```

4
doc/index.md
View File

@ -38,7 +38,6 @@ examples
Interfaces to calculators (VASP, QE, CRYSTAL, Abinit, TURBOMOLE) <interfaces>
command-options
input-output-files
hdf5_howto
auxiliary-tools
direct-solution
wigner-solution
@ -49,7 +48,8 @@ cutoff-pair
external-tools
phono3py-api
phono3py-load
tips
hdf5_howto
grid
citation
reference
changelog

479
doc/input-output-files.md
View File

@ -2,34 +2,30 @@
# Input / Output files
```{contents}
:depth: 3
:local:
```
The calculation results are written into files. Mostly the data are stored in
HDF5 format, therefore how to read the data from HDF5 files is also shown.
## Intermediate text files
The following files are not compatible with phonopy. But phonopy's `FORCE_SETS`
file can be created using phono3py command options from the following files. See
the detail at {ref}`file_format_compatibility`.
### `phono3py_disp.yaml`
This is created with `-d` option. See {ref}`create_displacements_option`.
## `phono3py_disp.yaml`
This is created with {ref}`-d <create_displacements_option>` or
{ref}`--rd <random_displacements_option>` option.
This file contains displacement dataset and crystal structure information.
Parameters for non-analytical term correction can be also included.
(input-output_files_FORCES_FC3)=
## `phono3py_params.yaml`
### `FORCES_FC3`
This is created with `--cf3` option. See {ref}`cf3_option`.
This is created with the combination of {ref}`--cf3 <cf3_option>` and {ref}`--sp
<sp_option>` options. This file contains displacement-force dataset and crystal
structure information. In addition, energies of supercells may be included in
the dataset. Parameters for non-analytical term correction can be also included.
There are two formats of `FORCES_FC3`. The type-I format is like that shown
below
(iofile_FORCES_FC3)=
## `FORCES_FC3`
This is created with {ref}`--cf3 <cf3_option>` option . There are two formats of
`FORCES_FC3`. The type-I format is like that shown below
```
# File: 1
@ -79,22 +75,13 @@ The type-II format is the same as
[phonopy's type-II format](https://phonopy.github.io/phonopy/input-files.html#type-2)
of `FORCE_SETS`.
### `FORCES_FC2`
## `FORCES_FC2`
This is created with `--cf2` option. See {ref}`cf2_option` and
{ref}`dim_fc2_option`.
This is created with {ref}`--cf2 <dim_fc2_option>` option. The file formats
(type-I and type-II) are same as those of `FORCES_FC3`.
The file formats (type-I and type-II) are same as those of `FORCES_FC3`.
## HDF5 files
### `kappa-*.hdf5`
See the detail at {ref}`kappa_hdf5_file`.
(fc3_hdf5_file)=
### `fc3.hdf5`
(iofile_fc3_hdf5)=
## `fc3.hdf5`
Third order force constants (in real space) are stored in
$\mathrm{eV}/\text{Angstrom}^3$.
@ -145,9 +132,8 @@ $$
So what you have to set is `--pa="0 1/4 1/4 1/4 0 1/4 1/4 1/4 0"`.
(fc2_hdf5_file)=
### `fc2.hdf5`
(iofile_fc2_hdf5)=
## `fc2.hdf5`
Second order force constants are stored in $\mathrm{eV}/\text{Angstrom}^2$.
@ -159,29 +145,434 @@ in the shape of:
```
against $\Phi_{\alpha\beta}(l\kappa, l'\kappa')$. More detail is similar to the
case for {ref}`fc3_hdf5_file`.
case for {ref}`iofile_fc3_hdf5`.
### `gamma-*.hdf5`
(iofile_kappa_hdf5)=
## `kappa-*.hdf5`
Files name, e.g. `kappa-m323220.hdf5`, is determined by some
specific options. `mxxx`, show the numbers of sampling
mesh. `sxxx` and `gxxx` appear optionally. `sxxx` gives the
smearing width in the smearing method for Brillouin zone integration
for phonon lifetime, and `gxxx` denotes the grid number. Using the
command option of `-o`, the file name can be modified slightly. For
example `-o nac` gives `kappa-m323220.nac.hdf5` to
memorize the option `--nac` was used.
### `mesh`
(Versions 1.10.11 or later)
The numbers of mesh points for reciprocal space sampling along
reciprocal axes, $a^*, b^*, c^*$.
### `frequency`
Phonon frequencies. The physical unit is THz, where THz
is in the ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band).
(iofile_kappa_hdf5_gamma)=
### `gamma`
Imaginary part of self energy of phonon bubble diagram (phonon-phonon
scattering). The physical unit is THz, where THz is in the ordinal frequency not
the angular frequency.
The array shape for all grid-points (irreducible q-points) is
(temperature, irreducible q-point, phonon band).
The array shape for a specific grid-point is
(temperature, phonon band).
Phonon lifetime ($\tau_\lambda=1/2\Gamma_\lambda(\omega_\lambda)$) may
be estimated from `gamma`. $2\pi$ has to be multiplied with
`gamma` values in the hdf5 file to convert the unit of ordinal
frequency to angular frequency. Zeros in `gamma` values mean that
those elements were not calculated such as for three acoustic modes at
$\Gamma$ point. The below is the copy-and-paste from the
previous section to show how to obtain phonon lifetime in pico
second:
```python
In [8]: g = f['gamma'][30]
In [9]: import numpy as np
In [10]: g = np.where(g > 0, g, -1)
In [11]: lifetime = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0)
```
### `gamma_isotope`
Isotope scattering of $1/2\tau^\mathrm{iso}_\lambda$.
The physical unit is same as that of gamma.
The array shape is same as that of frequency.
### `group_velocity`
Phonon group velocity, $\nabla_\mathbf{q}\omega_\lambda$. The
physical unit is $\text{THz}\cdot\text{Angstrom}$, where THz
is in the ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band, 3 = Cartesian coordinates).
### `heat_capacity`
Mode-heat-capacity defined by
$$
C_\lambda = k_\mathrm{B}
\left(\frac{\hbar\omega_\lambda}{k_\mathrm{B} T} \right)^2
\frac{\exp(\hbar\omega_\lambda/k_\mathrm{B}
T)}{[\exp(\hbar\omega_\lambda/k_\mathrm{B} T)-1]^2}.
$$
The physical unit is eV/K.
The array shape is (temperature, irreducible q-point, phonon band).
(iofile_kappa_hdf5_kappa)=
### `kappa`
Thermal conductivity tensor. The physical unit is W/m-K.
The array shape is (temperature, 6 = (xx, yy, zz, yz, xz, xy)).
### `mode-kappa`
Thermal conductivity tensors at k-stars (${}^*\mathbf{k}$):
$$
\sum_{\mathbf{q} \in {}^*\mathbf{k}} \kappa_{\mathbf{q}j}.
$$
The sum of this over ${}^*\mathbf{k}$ corresponding to
irreducible q-points divided by number of grid points gives
$\kappa$ ({ref}`iofile_kappa_hdf5_kappa`), e.g.,:
```python
kappa_xx_at_index_30 = mode_kappa[30, :, :, 0].sum()/ weight.sum()
```
Be careful that until version 1.12.7, mode-kappa values were divided
by number of grid points.
The physical unit is W/m-K. Each tensor element is the sum of tensor
elements on the members of ${}^*\mathbf{k}$, i.e., symmetrically
equivalent q-points by crystallographic point group and time reversal
symmetry.
The array shape is (temperature, irreducible q-point, phonon band, 6 =
(xx, yy, zz, yz, xz, xy)).
### `gv_by_gv`
Outer products of group velocities for k-stars
(${}^*\mathbf{k}$) for each irreducible q-point and phonon band
($j$):
$$
\sum_{\mathbf{q} \in {}^*\mathbf{k}} \mathbf{v}_{\mathbf{q}j} \otimes
\mathbf{v}_{\mathbf{q}j}.
$$
The physical unit is
$\text{THz}^2\cdot\text{Angstrom}^2$, where THz is in the
ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band, 6 = (xx, yy, zz,
yz, xz, xy)).
### `q-point`
Irreducible q-points in reduced coordinates.
The array shape is (irreducible q-point, 3 = reduced
coordinates in reciprocal space).
### `temperature`
Temperatures where thermal conductivities are calculated. The physical
unit is K.
### `weight`
Weights corresponding to irreducible q-points. Sum of weights equals to
the number of mesh grid points.
### `ave_pp`
Averaged phonon-phonon interaction $P_{\mathbf{q}j}$ in $\text{eV}^2$:
$$
P_{\mathbf{q}j} = \frac{1}{(3n_\mathrm{a})^2} \sum_{\lambda'\lambda''}
|\Phi_{\lambda\lambda'\lambda''}|^2.
$$
This is not going to be calculated in the RTA thermal coductivity
calculation mode by default. To calculate this, `--full-pp` option
has to be specified (see {ref}`full_pp_option`).
### `boundary_mfp`
A value specified by {ref}`boundary_mfp_option`. The physical unit is
micrometer.
When `--boundary-mfp` option is explicitly specified, its value is stored here.
### `kappa_unit_conversion`
This is used to convert the physical unit of lattice thermal
conductivity made of `heat_capacity`, `group_velocity`, and
`gamma`, to W/m-K. In the single mode relaxation time (SMRT) method,
a mode contribution to the lattice thermal conductivity is given by
$$
\kappa_\lambda = \frac{1}{V_0} C_\lambda \mathbf{v}_\lambda \otimes
\mathbf{v}_\lambda \tau_\lambda^{\mathrm{SMRT}}.
$$
For example, $\kappa_{\lambda,{xx}}$ is calculated by:
```python
In [1]: import h5py
In [2]: f = h5py.File("kappa-m111111.hdf5")
In [3]: kappa_unit_conversion = f['kappa_unit_conversion'][()]
In [4]: weight = f['weight'][:]
In [5]: heat_capacity = f['heat_capacity'][:]
In [6]: gv_by_gv = f['gv_by_gv'][:]
In [7]: gamma = f['gamma'][:]
In [8]: kappa_unit_conversion * heat_capacity[30, 2, 0] * gv_by_gv[2, 0] / (2 * gamma[30, 2, 0])
Out[8]:
array([ 1.02050241e+03, 1.02050241e+03, 1.02050241e+03,
4.40486382e-15, 0.00000000e+00, -4.40486382e-15])
In [9]: f['mode_kappa'][30, 2, 0]
Out[9]:
array([ 1.02050201e+03, 1.02050201e+03, 1.02050201e+03,
4.40486209e-15, 0.00000000e+00, -4.40486209e-15])
```
(iofile_kappa_hdf5_gamma_NU)=
### `gamma_N` and `gamma_U`
The data are stored in `kappa-mxxx(-gx-sx-sdx).hdf5` file and accessed by
`gamma_N` and `gamma_U` keys. The shape of the arrays is the same as that of
`gamma` (see {ref}`iofile_kappa_hdf5_gamma`). An example (Si-PBEsol) is shown
below:
```bash
% phono3py-load --mesh 11 11 11 --fc3 --fc2 --br --nu
...
% ipython
```
```python
In [1]: import h5py
In [2]: f = h5py.File("kappa-m111111.hdf5", 'r')
In [3]: list(f)
Out[3]:
['frequency',
'gamma',
'gamma_N',
'gamma_U',
'group_velocity',
'gv_by_gv',
'heat_capacity',
'kappa',
'kappa_unit_conversion',
'mesh',
'mode_kappa',
'qpoint',
'temperature',
'weight']
In [4]: f['gamma'].shape
Out[4]: (101, 56, 6)
In [5]: f['gamma_N'].shape
Out[5]: (101, 56, 6)
In [6]: f['gamma_U'].shape
Out[6]: (101, 56, 6)
```
## `gamma-*.hdf5`
Imaginary parts of self energies at harmonic phonon frequencies
($\Gamma_\lambda(\omega_\lambda)$ = half linewidths) are stored in THz. See
{ref}`write_gamma_option`.
### `gamma_detail-*.hdf5`
(iofile_gamma_detail_hdf5)=
## `gamma_detail-*.hdf5`
Q-point triplet contributions to imaginary parts of self energies at phonon
frequencies (half linewidths) are stored in THz. See
{ref}`write_detailed_gamma_option`.
{ref}`--write-gamma-detail <write_detailed_gamma_option>` option.
## Simple text file
In the output file in hdf5, following keys are used to extract the detailed
information.
### `gammas-*.dat`
```{table}
| dataset | Array shape |
| --------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------ |
| gamma_detail for `--ise` | (temperature, sampling frequency point, symmetry reduced set of triplets at given grid point, band1, band2, band3) in THz (without $2\pi$) |
| gamma_detail for `--br` | (temperature, symmetry reduced set of triplets at gvien grid point, band1, band2, band3) in THz (without $2\pi$) |
| mesh | Numbers of sampling mesh along reciprocal axes. |
| frequency_point for `--ise` | Sampling frequency points in THz (without $2\pi$), i.e., $\omega$ in $\Gamma_\lambda(\omega)$ |
| temperature | (temperature,), Temperatures in K |
| triplet | (symmetry reduced set of triplets at given grid point, 3), Triplets are given by the grid point indices (see below). |
| weight | (symmetry reduced set of triplets at given grid point,), Weight of each triplet to imaginary part of self energy |
| triplet_all | (triplets at given grid point, 3), symmetry non-reduced version of the triplet information. |
```
See {ref}`grid_triplets` to recover the q-points of each triplet.
Imaginary part of self energy (linewidth/2) is recovered by the following
script:
```python
import h5py
import numpy as np
gd = h5py.File("gamma_detail-mxxx-gx.hdf5")
temp_index = 30 # index of temperature
temperature = gd['temperature'][temp_index]
gamma_tp = gd['gamma_detail'][:].sum(axis=-1).sum(axis=-1)
weight = gd['weight'][:]
gamma = np.dot(weight, gamma_tp[temp_index])
```
For example, for `--br`, this `gamma` gives $\Gamma_\lambda(\omega_\lambda)$ of
the band indices at the grid point indicated by $\lambda$ at the temperature of
index 30. If any bands are degenerated, those `gamma` in
`kappa-mxxx-gx(-sx-sdx).hdf5` or `gamma-mxxx-gx(-sx-sdx).hdf5` type file are
averaged, but the `gamma` obtained here in this way are not symmetrized. Apart
from this symmetrization, the values must be equivalent between them.
To understand each contribution of triptle to imaginary part of self energy,
reading `phonon-mxxx.hdf5` is useful (see {ref}`write_phonon_option`). For
example, phonon triplets of three phonon scatterings are obtained by
```python
import h5py
import numpy as np
gd = h5py.File("gamma_detail-mxxx-gx.hdf5", 'r')
ph = h5py.File("phonon-mxxx.hdf5", 'r')
gp1 = gd['grid_point'][()]
triplets = gd['triplet'][:] # Sets of (gp1, gp2, gp3) where gp1 is fixed
mesh = gd['mesh'][:]
grid_address = ph['grid_address'][:]
q_triplets = grid_address[triplets] / mesh.astype('double') # For conventional regular grid
# Phonons of triplets[2]
phonon_tp = [(ph['frequency'][i], ph['eigenvector'][i]) for i in triplets[2]]
# Fractions of contributions of triplets at this grid point and temperature index 30
gamma_sum_over_bands = np.dot(weight, gd['gamma_detail'][30].sum(axis=-1).sum(axis=-1).sum(axis=-1))
contrib_tp = [gd['gamma_detail'][30, i].sum() / gamma_sum_over_bands for i in range(len(weight))]
np.dot(weight, contrib_tp) # is one
```
(iofile_phonon_hdf5)=
## `phonon-*.hdf5`
Contents of `phonon-mxxx.hdf5` are watched by:
```python
In [1]: import h5py
In [2]: f = h5py.File("phonon-m111111.hdf5")
In [3]: list(f)
Out[3]:
['eigenvector',
'frequency',
'grid_address',
'ir_grid_points',
'ir_grid_weights',
'mesh',
'version']
In [4]: f['mesh'][:]
Out[4]: array([11, 11, 11])
In [5]: f['grid_address'].shape
Out[5]: (1367, 3)
In [6]: f['frequency'].shape
Out[6]: (1367, 6)
In [7]: f['eigenvector'].shape
Out[7]: (1367, 6, 6)
In [8]: f['ir_grid_points'].shape
Out[8]: (56,)
```
The first axis of `ph['grid_address']`, `ph['frequency']`, and
`ph['eigenvector']` corresponds to the number of q-points where phonons are
calculated. Here the number of phonons may not be equal to product of mesh
numbers ($1367 \neq 11^3$). This is because all q-points on Brillouin zone
boundary are included, i.e., even if multiple q-points are translationally
equivalent, those phonons are stored separately though these phonons are
physically equivalent within the equations employed in phono3py. Here Brillouin
zone is defined by Wigner–Seitz cell of reciprocal primitive basis vectors. This
is convenient to categorize phonon triplets into Umklapp and Normal scatterings
based on the Brillouin zone.
## `pp-*.hdf5`
This file contains phonon-phonon interaction strength
$\bigl|\Phi_{\lambda\lambda'\lambda''}\bigl|^2$. To use the data in this
file, it is recommended to generate with `--full-pp` option because the data
structure to access becomes simpler.
```bash
% phono3py-load phono3py.yaml --gp 5 --br --mesh 11 11 11 --write-pp --full-pp
```
```python
In [1]: import h5py
In [2]: f = h5py.File("pp-m111111-g5.hdf5")
In [3]: list(f)
Out[3]: ['pp', 'triplet', 'triplet_all', 'version', 'weight']
In [4]: f['pp'].shape
Out[4]: (146, 6, 6, 6)
```
Indices of the `pp` array are (symmetry reduced set of triplets at given grid
point, band1, band2, band3), and the values are given in $\text{eV}^2$. See
{ref}`grid_triplets` to recover the q-points of each triplet.
Except for `pp`, all the other information are equivalent to those found in
{ref}`iofile_gamma_detail_hdf5`.
## `gammas-*.dat`
Imaginary parts of self energies with respect to frequency
$\Gamma_\lambda(\omega)$ are stored in THz. See {ref}`ise_option`.
### `jdos-*.dat`
## `jdos-*.dat`
Joint densities of states are stored in Thz. See {ref}`jdos_option`.
### `linewidth-*.dat`
## `linewidth-*.dat`

2
doc/phono3py-api.md
View File

@ -148,7 +148,7 @@ shape of `(num_supercells, num_atoms_in_supercell, 3)`. In the above example,
the array shape is `(1254, 72, 3)`.
If the calculated force sets are stored in the
{ref}`input-output_files_FORCES_FC3` file, the numpy array of `forces` is
{ref}`iofile_FORCES_FC3` file, the numpy array of `forces` is
obtained by
```python

113
doc/tips.md
View File

@ -1,113 +0,0 @@
(tips)=
# Tips
```{contents}
:depth: 2
:local:
```
(brillouinzone_sum)=
## Brillouin zone summation
Brillouin zone (BZ) summations appear at different two points in
phonon lifetime calculation. First it is used for the Fourier
transform of force constants, and second to obtain imaginary part of
phonon-self-energy. For the summation, usually uniform sampling meshes
are employed. To obtain more accurate result, it is always better to
use denser meshes. But the denser mesh requires more computational
demand.
The second BZ summation contains delta functions. In
phono3py calculation, a linear tetrahedron method ({ref}`thm option <thm_option>`, default option) and a smearing method ({ref}`sigma option <sigma_option>`) can be used for this BZ
integration. In most cases, the tetrahedron method is better. Especially in high
thermal conductivity materials, the smearing method results in
underestimation of lattice thermal conductivity.
The figure below shows Si thermal conductivity convergence with
respect to number of mesh points along an axis from n=19 to 65. This
is calculated with RTA and the linear tetrahedron method. Within the
methods and phono3py implementation, it is converging at around n=55,
however this computational demand is not trivial. As far as observing
this result, an extrapolation to $1/n \rightarrow 0$ seems not a
good idea, since it gives overestimation in the case of this Si
example. This plot tells that we have to decide how much value is
acceptable as lattice thermal conductivity value. Therefore it
important to describe the number of sampling mesh and method of BZ
integration to let other people reproduce the computational results.
```{image} Si-convergence.png
:width: 25%
```
In case the smearing method is necessary to use, the convergence of
q-point mesh together with smearing width has to be checked
carefully. Since smearing parameter is used to approximate delta
functions, small `sigma` value is better to describe the detailed
structure of three-phonon-space, however it requires a denser sampling
mesh to converge the result. To check the convergence with respect to
the `sigma` value, multiple sigma values can be set. This can be
computationally efficient, since it is avoided to re-calculate
phonon-phonon interaction strength for different `sigma` values in
this case. Convergence with respect to the sampling mesh and smearing
parameter strongly depends on materials. For Si example, a
$20\times 20\times 20$ sampling mesh (or 8000 reducible sampling
points) and 0.1 THz smearing value for reciprocal of the volume of an
atom may be a good starting choice. The tetrahedron method requires no
such parameter as the smearing width.
## Importance of numerical quality of force constants
Third-order force constants (fc3) are much weaker to numerical noise
of a force calculator than second-order force constants
(fc2). Therefore supercell force calculations have to be done very
carefully.
Numerical quality of forces given by force calculators is the most
important factor for the numerical quality of lattice thermal
conductivity calculation. We may be able to apply symmetry constraints
to force constants, however even if force constants fulfill those
symmetries, the numerical quality of force constants is not guaranteed
since elements of force constants just suffice the symmetries but most
of those intensities are not constrained.
It is important to use the best possible force calculator in the
possibly best way. The knowledge of the force calculator from the
theory and method to the practical usage is required to obtain
good results of lattice thermal conductivity calculation.
In the following, a few things that may be good to know are
presented.
### A practical way to check lattice thermal conductivity result
Some feeling whether our calculation result is OK or not may be
obtained by comparing lattice thermal conductivities calculated with
and without {ref}`symmetrizations of force constants <symmetrization_option>`. If they are enough different, e.g., more
than twice different, it is better to re-consider about the force
calculation. In the case of DFT calculations, the choice of input
settings such as k-point sampling mesh, plane-wave energy cutoff, and
exchange-correlational potential, etc, should be reconsidered.
### Displacement distance of atoms
The phono3py default displacement distance is 0.03
$\text{Angstrom}$. In some cases, accurate result may not be obtained
due to the numerical noise of the force calculator. Usually increasing
the displacement distance by the {ref}`amplitude option <amplitude_option>` reduces the numerical noise, but as its drawback
higher order anharmonicity is involved (renormalized) into fc3 and fc2.
(file_format_compatibility)=
## File format compatibility with phonopy
- `FORCES_FC3` and `FORCES_FC2` are not
compatible with phonopy's `FORCE_SETS`.
- `FORCE_SETS` can be created using {ref}`--cfs <cfs_option>` from
`FORCES_FC3` and `phono3py_disp.yaml` or `FORCES_FC2` and
`phono3py_disp.yaml` (needs to specify `--dim-fc2`).
- `FORCES_FC2` can be created using {ref}`--fs2f2 <fs2f2_option>` from `FORCE_SETS`.
- `fc2.hdf5` can be used in phonopy in the `hdf5` mode when it is
renamed to `force_constants.hdf5`. In the previous combinations of
phonopy and phono3py, depending on the physical unit of force
constants of calculators, the direct compatibility is not guaranteed.

135
phono3py/api_phono3py.py
View File

@ -38,9 +38,10 @@ from __future__ import annotations
import copy
import warnings
from collections.abc import Sequence
from typing import Literal, Optional, Union
from typing import Literal, Optional, Union, cast
import numpy as np
from numpy.typing import NDArray
from phonopy.harmonic.displacement import (
directions_to_displacement_dataset,
get_least_displacements,
@ -58,6 +59,7 @@ from phonopy.interface.mlp import PhonopyMLP
from phonopy.interface.pypolymlp import (
PypolymlpParams,
)
from phonopy.interface.symfc import SymfcFCSolver
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import (
@ -73,7 +75,7 @@ from phonopy.structure.symmetry import Symmetry
from phono3py.conductivity.init_direct_solution import get_thermal_conductivity_LBTE
from phono3py.conductivity.init_rta import get_thermal_conductivity_RTA
from phono3py.interface.fc_calculator import get_fc3_solver
from phono3py.interface.fc_calculator import FC3Solver
from phono3py.interface.phono3py_yaml import Phono3pyYaml
from phono3py.phonon.grid import BZGrid
from phono3py.phonon3.dataset import forces_in_dataset
@ -249,10 +251,10 @@ class Phono3py:
)
else:
self._phonon_supercell_matrix = None
self._supercell = None
self._primitive = None
self._phonon_supercell = None
self._phonon_primitive = None
self._supercell: Supercell
self._primitive: Primitive
self._phonon_supercell: Supercell
self._phonon_primitive: Primitive
self._build_supercell()
self._build_primitive_cell()
self._build_phonon_supercell()
@ -267,9 +269,9 @@ class Phono3py:
self._bz_grid = None
# Set supercell, primitive, and phonon supercell symmetries
self._symmetry = None
self._primitive_symmetry = None
self._phonon_supercell_symmetry = None
self._symmetry: Symmetry
self._primitive_symmetry: Symmetry
self._phonon_supercell_symmetry: Symmetry
self._search_symmetry()
self._search_primitive_symmetry()
self._search_phonon_supercell_symmetry()
@ -298,6 +300,8 @@ class Phono3py:
# Force constants
self._fc2 = None
self._fc3 = None
self._fc3_nonzero_indices = None # available only symfc
self._fc3_cutoff = None # available only symfc
# MLP
self._mlp = None
@ -332,10 +336,10 @@ class Phono3py:
return self._calculator
@property
def fc3(self) -> Optional[np.ndarray]:
def fc3(self) -> NDArray | None:
"""Setter and getter of third order force constants (fc3).
ndarray
ndarray, optional
fc3 shape is either (supercell, supercell, supercell, 3, 3, 3) or
(primitive, supercell, supercell, 3, 3, 3),
where 'supercell' and 'primitive' indicate number of atoms in
@ -349,7 +353,26 @@ class Phono3py:
self._fc3 = fc3
@property
def fc2(self) -> Optional[np.ndarray]:
def fc3_nonzero_indices(self) -> NDArray | None:
"""Setter and getter of non-zero indices of fc3.
ndarray, optional
Non-zero indices of fc3.
"""
return self._fc3_nonzero_indices
@fc3_nonzero_indices.setter
def fc3_nonzero_indices(self, fc3_nonzero_indices):
self._fc3_nonzero_indices = fc3_nonzero_indices
@property
def fc3_cutoff(self) -> float | None:
"""Return cutoff value of fc3."""
return self._fc3_cutoff
@property
def fc2(self) -> NDArray | None:
"""Setter and getter of second order force constants (fc2).
ndarray
@ -366,7 +389,7 @@ class Phono3py:
self._fc2 = fc2
@property
def force_constants(self) -> Optional[np.ndarray]:
def force_constants(self) -> NDArray | None:
"""Return fc2. This is same as the getter attribute `fc2`."""
return self.fc2
@ -537,7 +560,7 @@ class Phono3py:
return self._phonon_supercell_symmetry
@property
def supercell_matrix(self) -> np.ndarray:
def supercell_matrix(self) -> NDArray:
"""Return transformation matrix to supercell cell from unit cell.
ndarray
@ -548,7 +571,7 @@ class Phono3py:
return self._supercell_matrix
@property
def phonon_supercell_matrix(self) -> Optional[np.ndarray]:
def phonon_supercell_matrix(self) -> NDArray | None:
"""Return transformation matrix to phonon supercell from unit cell.
ndarray
@ -559,7 +582,7 @@ class Phono3py:
return self._phonon_supercell_matrix
@property
def primitive_matrix(self) -> np.ndarray:
def primitive_matrix(self) -> NDArray:
"""Return transformation matrix to primitive cell from unit cell.
ndarray
@ -730,10 +753,10 @@ class Phono3py:
self._mlp = mlp
@property
def band_indices(self) -> list[np.ndarray]:
def band_indices(self) -> list[NDArray]:
"""Setter and getter of band indices.
list[np.ndarray]
list[NDArray]
List of band indices specified to select specific bands
to computer ph-ph interaction related properties.
@ -753,7 +776,7 @@ class Phono3py:
self._band_indices_flatten = np.hstack(self._band_indices).astype("int64")
@property
def masses(self) -> np.ndarray:
def masses(self) -> NDArray:
"""Setter and getter of atomic masses of primitive cell."""
return self._primitive.masses
@ -814,7 +837,7 @@ class Phono3py:
return self._bz_grid.D_diag
@mesh_numbers.setter
def mesh_numbers(self, mesh_numbers: Union[int, float, Sequence, np.ndarray]):
def mesh_numbers(self, mesh_numbers: Union[int, float, Sequence, NDArray]):
self._set_mesh_numbers(mesh_numbers)
@property
@ -1121,7 +1144,7 @@ class Phono3py:
`None` will choose one of them automatically.
"""
if self.mesh_numbers is None:
if self._bz_grid is None:
msg = "Phono3py.mesh_numbers of instance has to be set."
raise RuntimeError(msg)
@ -1134,6 +1157,7 @@ class Phono3py:
self._bz_grid,
self._primitive_symmetry,
fc3=self._fc3,
fc3_nonzero_indices=self._fc3_nonzero_indices,
band_indices=self._band_indices_flatten,
constant_averaged_interaction=constant_averaged_interaction,
frequency_factor_to_THz=self._frequency_factor_to_THz,
@ -1318,15 +1342,15 @@ class Phono3py:
number_of_snapshots == "auto" or number_of_snapshots > 0
):
if number_of_snapshots == "auto":
from phonopy.interface.symfc import SymfcFCSolver
if cutoff_pair_distance is None:
options = None
else:
options = {"cutoff": {3: cutoff_pair_distance}}
_number_of_snapshots = (
SymfcFCSolver(
self._supercell, self._symmetry, options=options
self._supercell,
symmetry=self._symmetry,
options=options,
).estimate_numbers_of_supercells(orders=[3])[3]
* 2
)
@ -1340,6 +1364,8 @@ class Phono3py:
random_seed=random_seed,
max_distance=max_distance,
)
if cutoff_pair_distance is not None:
self._dataset["cutoff_distance"] = cutoff_pair_distance
else:
direction_dataset = get_third_order_displacements(
self._supercell,
@ -1426,11 +1452,9 @@ class Phono3py:
number_of_snapshots == "auto" or number_of_snapshots > 0
):
if number_of_snapshots == "auto":
from phonopy.interface.symfc import SymfcFCSolver
_number_of_snapshots = (
SymfcFCSolver(
self._supercell, self._symmetry
self._supercell, symmetry=self._symmetry
).estimate_numbers_of_supercells(orders=[2])[2]
* 2
)
@ -1486,14 +1510,16 @@ class Phono3py:
Options for external force constants calculator.
"""
fc_solver = get_fc3_solver(
fc_solver_name = fc_calculator if fc_calculator is not None else "traditional"
fc_solver = FC3Solver(
fc_solver_name,
self._supercell,
self._primitive,
self._dataset,
fc_calculator=fc_calculator,
fc_calculator_options=fc_calculator_options,
is_compact_fc=is_compact_fc,
symmetry=self._symmetry,
dataset=self._dataset,
is_compact_fc=is_compact_fc,
primitive=self._primitive,
orders=[2, 3],
options=fc_calculator_options,
log_level=self._log_level,
)
fc2 = fc_solver.force_constants[2]
@ -1512,8 +1538,25 @@ class Phono3py:
if self._fc2 is None:
symmetrize_force_constants(fc2)
# Set fc2 and fc3
self._fc3 = fc3
self._fc3_nonzero_indices = None
if fc_calculator == "symfc":
symfc_solver = cast(SymfcFCSolver, fc_solver.fc_solver)
fc3_nonzero_elems = symfc_solver.get_nonzero_atomic_indices_fc3()
options = symfc_solver.options
if options is not None and "cutoff" in options:
self._fc3_cutoff = options["cutoff"].get(3, None)
if fc3_nonzero_elems is not None:
if is_compact_fc:
self._fc3_nonzero_indices = np.array(
fc3_nonzero_elems[self._primitive.p2s_map],
dtype="byte",
order="C",
)
else:
self._fc3_nonzero_indices = np.array(
fc3_nonzero_elems, dtype="byte", order="C"
)
# fc2 as obtained above will not be set when "|" in fc-calculator setting.
if fc_calculator is not None and "|" in fc_calculator:
@ -2349,7 +2392,9 @@ class Phono3py:
# private methods #
###################
def _search_symmetry(self):
self._symmetry = Symmetry(self._supercell, self._symprec, self._is_symmetry)
self._symmetry = Symmetry(
self._supercell, symprec=self._symprec, is_symmetry=self._is_symmetry
)
def _search_primitive_symmetry(self):
self._primitive_symmetry = Symmetry(
@ -2357,7 +2402,7 @@ class Phono3py:
)
if len(self._symmetry.pointgroup_operations) != len(
self._primitive_symmetry.pointgroup_operations
): # noqa E129
):
print(
"Warning: point group symmetries of supercell and primitive"
"cell are different."
@ -2368,12 +2413,14 @@ class Phono3py:
self._phonon_supercell_symmetry = self._symmetry
else:
self._phonon_supercell_symmetry = Symmetry(
self._phonon_supercell, self._symprec, self._is_symmetry
self._phonon_supercell,
symprec=self._symprec,
is_symmetry=self._is_symmetry,
)
def _build_supercell(self):
self._supercell = get_supercell(
self._unitcell, self._supercell_matrix, self._symprec
self._unitcell, self._supercell_matrix, symprec=self._symprec
)
def _build_primitive_cell(self):
@ -2405,7 +2452,7 @@ class Phono3py:
self._phonon_supercell = self._supercell
else:
self._phonon_supercell = get_supercell(
self._unitcell, self._phonon_supercell_matrix, self._symprec
self._unitcell, self._phonon_supercell_matrix, symprec=self._symprec
)
def _build_phonon_primitive_cell(self):
@ -2498,7 +2545,9 @@ class Phono3py:
self._supercells_with_displacements = supercells
def _get_primitive_cell(self, supercell, supercell_matrix, primitive_matrix):
def _get_primitive_cell(
self, supercell, supercell_matrix, primitive_matrix
) -> Primitive:
inv_supercell_matrix = np.linalg.inv(supercell_matrix)
if primitive_matrix is None:
t_mat = inv_supercell_matrix
@ -2516,7 +2565,7 @@ class Phono3py:
def _set_mesh_numbers(
self,
mesh: Union[int, float, Sequence, np.ndarray],
mesh: Union[int, float, Sequence, NDArray],
):
# initialization related to mesh
self._interaction = None
@ -2562,7 +2611,7 @@ class Phono3py:
def _get_forces_energies(
self, target: Literal["forces", "supercell_energies"]
) -> Optional[np.ndarray]:
) -> NDArray | None:
"""Return fc3 forces and supercell energies.
Return None if tagert data is not found rather than raising exception.
@ -2627,7 +2676,7 @@ class Phono3py:
def _get_phonon_forces_energies(
self, target: Literal["forces", "supercell_energies"]
) -> Optional[np.ndarray]:
) -> NDArray | None:
"""Return fc2 forces and supercell energies.
Return None if tagert data is not found rather than raising exception.

709
phono3py/conductivity/base.py
View File

@ -34,18 +34,24 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import textwrap
import warnings
from abc import ABC, abstractmethod
from typing import List, Optional, Tuple, Union
from typing import Optional
import numpy as np
from numpy.typing import NDArray
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.phonon.thermal_properties import mode_cv
from phonopy.physical_units import get_physical_units
from phono3py.other.isotope import Isotope
from phono3py.phonon.grid import get_grid_points_by_rotations, get_ir_grid_points
from phono3py.phonon.grid import (
get_grid_points_by_rotations,
get_ir_grid_points,
get_qpoints_from_bz_grid_points,
)
from phono3py.phonon3.collision_matrix import CollisionMatrix
from phono3py.phonon3.imag_self_energy import ImagSelfEnergy
from phono3py.phonon3.interaction import Interaction
@ -63,12 +69,135 @@ def get_unit_to_WmK() -> float:
return unit_to_WmK
class HeatCapacityMixIn:
"""Heat capacity mix-in.
def get_multiplicity_at_q(
gp: int,
pp: Interaction,
point_operations: np.ndarray,
) -> int:
"""Return multiplicity (order of site-symmetry) of q-point."""
q = get_qpoints_from_bz_grid_points(gp, pp.bz_grid)
reclat = np.linalg.inv(pp.primitive.cell)
multi = 0
for q_rot in [np.dot(r, q) for r in point_operations]:
diff = q - q_rot
diff -= np.rint(diff)
dist = np.linalg.norm(np.dot(reclat, diff))
if dist < pp.primitive_symmetry.tolerance:
multi += 1
return multi
Used by other mix-in.
def get_kstar_order(
grid_weight: int, multi: int, point_operations: np.ndarray, verbose: bool = False
) -> int:
"""Return order (number of arms) of kstar.
multi : int
Multiplicity of grid point.
"""
order_kstar = len(point_operations) // multi
if order_kstar != grid_weight:
if verbose:
text = (
"Number of elements in k* is unequal "
"to number of equivalent grid-points. "
"This means that the mesh sampling grids break "
"symmetry. Please check carefully "
"the convergence over grid point densities."
)
msg = textwrap.fill(
text, initial_indent=" ", subsequent_indent=" ", width=70
)
print("*" * 30 + "Warning" + "*" * 30)
print(msg)
print("*" * 67)
return order_kstar
def get_heat_capacities(
grid_point: int,
pp: Interaction,
temperatures: NDArray[np.float64],
):
"""Return mode heat capacity.
cv returned should be given to self._cv by
self._cv[:, i_data, :] = cv
"""
if not pp.phonon_all_done:
raise RuntimeError(
"Phonon calculation has not been done yet. "
"Run phono3py.run_phonon_solver() before this method."
)
frequencies, _, _ = pp.get_phonons()
freqs = (
frequencies[grid_point][pp.band_indices] # type: ignore
* get_physical_units().THzToEv
)
cutoff = pp.cutoff_frequency * get_physical_units().THzToEv
cv = np.zeros((len(temperatures), len(freqs)), dtype="double")
# x=freq/T has to be small enough to avoid overflow of exp(x).
# x < 100 is the hard-corded criterion.
# Otherwise just set 0.
for i, f in enumerate(freqs):
if f > cutoff:
condition = f < 100 * temperatures * get_physical_units().KB
cv[:, i] = np.where(
condition,
mode_cv(np.where(condition, temperatures, 10000), f), # type: ignore
0,
)
return cv
class ConductivityComponentsBase(ABC):
"""Base class of ConductivityComponents."""
def __init__(
self,
pp: Interaction,
grid_points: NDArray[np.int64],
grid_weights: NDArray[np.int64],
point_operations: NDArray[np.int64],
rotations_cartesian: NDArray[np.int64],
temperatures: Optional[NDArray[np.float64]] = None,
average_gv_over_kstar: bool = False,
is_kappa_star: bool = True,
gv_delta_q: float | None = None,
is_reducible_collision_matrix: bool = False,
log_level: int = 0,
):
"""Init method.
Parameters
----------
gv_delta_q : float, optional, default is None, # for group velocity
With non-analytical correction, group velocity is calculated
by central finite difference method. This value gives the distance
in both directions in 1/Angstrom. The default value will be 1e-5.
"""
self._pp = pp
self._grid_points = grid_points
self._grid_weights = grid_weights
self._point_operations = point_operations
self._rotations_cartesian = rotations_cartesian
self._temperatures = temperatures
self._average_gv_over_kstar = average_gv_over_kstar
self._gv_delta_q = gv_delta_q
self._is_kappa_star = is_kappa_star
self._is_reducible_collision_matrix = is_reducible_collision_matrix
self._log_level = log_level
self._gv: np.ndarray
self._cv: np.ndarray
self._num_sampling_grid_points = 0
@property
def mode_heat_capacities(self):
@ -79,87 +208,6 @@ class HeatCapacityMixIn:
"""
return self._cv
def get_mode_heat_capacities(self):
"""Return mode heat capacity at constant volume at grid points.
Grid points are those at mode kappa are calculated.
"""
warnings.warn(
"Use attribute, Conductivity.mode_heat_capacities "
"instead of Conductivity.get_mode_heat_capacities().",
DeprecationWarning,
stacklevel=2,
)
return self.mode_heat_capacities
def _set_cv(self, i_gp, i_data):
"""Set mode heat capacity.
The array has to be allocated somewhere out of the mix-in.
self._cv = np.zeros(
(num_temp, num_grid_points, num_band0), order="C", dtype="double"
)
"""
grid_point = self._grid_points[i_gp]
freqs = (
self._frequencies[grid_point][self._pp.band_indices]
* get_physical_units().THzToEv
)
cutoff = self._pp.cutoff_frequency * get_physical_units().THzToEv
cv = np.zeros((len(self._temperatures), len(freqs)), dtype="double")
# x=freq/T has to be small enough to avoid overflow of exp(x).
# x < 100 is the hard-corded criterion.
# Otherwise just set 0.
for i, f in enumerate(freqs):
if f > cutoff:
condition = f < 100 * self._temperatures * get_physical_units().KB
cv[:, i] = np.where(
condition,
mode_cv(np.where(condition, self._temperatures, 10000), f),
0,
)
self._cv[:, i_data, :] = cv
class ConductivityMixIn(HeatCapacityMixIn):
"""Thermal conductivity mix-in.
Used by ConductivityRTA and ConductivityLBTE.
"""
@property
def kappa(self):
"""Return kappa."""
return self._kappa
def get_kappa(self):
"""Return kappa."""
warnings.warn(
"Use attribute, Conductivity.kappa instead of Conductivity.get_kappa().",
DeprecationWarning,
stacklevel=2,
)
return self.kappa
@property
def mode_kappa(self):
"""Return mode_kappa."""
return self._mode_kappa
def get_mode_kappa(self):
"""Return mode_kappa."""
warnings.warn(
"Use attribute, Conductivity.mode_kappa "
"instead of Conductivity.get_mode_kappa().",
DeprecationWarning,
stacklevel=2,
)
return self.mode_kappa
@property
def group_velocities(self):
"""Return group velocities at grid points.
@ -169,36 +217,96 @@ class ConductivityMixIn(HeatCapacityMixIn):
"""
return self._gv
def get_group_velocities(self):
"""Return group velocities at grid points.
Grid points are those at mode kappa are calculated.
"""
warnings.warn(
"Use attribute, Conductivity.group_velocities "
"instead of Conductivity.get_group_velocities().",
DeprecationWarning,
stacklevel=2,
)
return self.group_velocities
@property
def gv_delta_q(self):
"""Return delta q for group velocity."""
return self._gv_delta_q
@property
def gv_by_gv(self):
"""Return gv_by_gv at grid points where mode kappa are calculated."""
return self._gv_sum2
def number_of_sampling_grid_points(self):
"""Return number of grid points.
def get_gv_by_gv(self):
"""Return gv_by_gv at grid points where mode kappa are calculated."""
warnings.warn(
"Use attribute, Conductivity.gv_by_gv "
"instead of Conductivity.get_gv_by_gv().",
DeprecationWarning,
stacklevel=2,
This is calculated by the sum of numbers of arms of k-start.
"""
return self._num_sampling_grid_points
def set_heat_capacities(self, i_gp, i_data):
"""Set heat capacity at grid point and at data location."""
if self._temperatures is None:
raise RuntimeError(
"Temperatures have not been set yet. "
"Set temperatures before this method."
)
cv = get_heat_capacities(self._grid_points[i_gp], self._pp, self._temperatures)
self._cv[:, i_data, :] = cv
@abstractmethod
def set_velocities(self, i_gp, i_data):
"""Set velocities at grid point and at data location."""
raise NotImplementedError()
def _allocate_values(self):
if self._temperatures is None:
raise RuntimeError(
"Temperatures have not been set yet. "
"Set temperatures before this method."
)
num_band0 = len(self._pp.band_indices)
if self._is_reducible_collision_matrix:
num_grid_points = np.prod(self._pp.mesh_numbers)
else:
num_grid_points = len(self._grid_points)
num_temp = len(self._temperatures)
self._cv = np.zeros(
(num_temp, num_grid_points, num_band0), order="C", dtype="double"
)
return self.gv_by_gv
self._gv = np.zeros((num_grid_points, num_band0, 3), order="C", dtype="double")
def _init_velocity(self, gv_delta_q):
class ConductivityComponents(ConductivityComponentsBase):
"""Thermal conductivity components.
Used by ConductivityRTA and ConductivityLBTE.
"""
def __init__(
self,
pp: Interaction,
grid_points: NDArray[np.int64],
grid_weights: NDArray[np.int64],
point_operations: NDArray[np.int64],
rotations_cartesian: NDArray[np.int64],
temperatures: Optional[NDArray[np.float64]] = None,
average_gv_over_kstar: bool = False,
is_kappa_star: bool = True,
gv_delta_q: Optional[float] = None,
is_reducible_collision_matrix: bool = False,
log_level: int = 0,
):
"""Init method."""
super().__init__(
pp,
grid_points,
grid_weights,
point_operations,
rotations_cartesian,
temperatures=temperatures,
average_gv_over_kstar=average_gv_over_kstar,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_reducible_collision_matrix=is_reducible_collision_matrix,
log_level=log_level,
)
self._gv_by_gv: np.ndarray
if self._pp.dynamical_matrix is None:
raise RuntimeError("Interaction.init_dynamical_matrix() has to be called.")
self._velocity_obj = GroupVelocity(
self._pp.dynamical_matrix,
q_length=gv_delta_q,
@ -206,7 +314,16 @@ class ConductivityMixIn(HeatCapacityMixIn):
frequency_factor_to_THz=self._pp.frequency_factor_to_THz,
)
def _set_velocities(self, i_gp, i_data):
if self._temperatures is not None:
self._allocate_values()
@property
def gv_by_gv(self):
"""Return gv_by_gv at grid points where mode kappa are calculated."""
return self._gv_by_gv
def set_velocities(self, i_gp, i_data):
"""Set group velocities at grid point and at data location."""
self._gv[i_data] = self._get_gv(i_gp)
self._set_gv_by_gv(i_gp, i_data)
@ -222,8 +339,11 @@ class ConductivityMixIn(HeatCapacityMixIn):
unique_gps = np.unique(gps_rotated)
gvs = {}
for bz_gp in unique_gps.tolist(): # To convert to int type.
self._velocity_obj.run([self._get_qpoint_from_gp_index(bz_gp)])
for bz_gp in unique_gps:
self._velocity_obj.run(
[get_qpoints_from_bz_grid_points(bz_gp, self._pp.bz_grid)]
)
assert self._velocity_obj.group_velocities is not None
gvs[bz_gp] = self._velocity_obj.group_velocities[
0, self._pp.band_indices, :
]
@ -232,7 +352,10 @@ class ConductivityMixIn(HeatCapacityMixIn):
gv += np.dot(gvs[bz_gp], r) # = dot(r_inv, gv)
return gv / len(self._point_operations)
else:
self._velocity_obj.run([self._get_qpoint_from_gp_index(irgp)])
self._velocity_obj.run(
[get_qpoints_from_bz_grid_points(irgp, self._pp.bz_grid)]
)
assert self._velocity_obj.group_velocities is not None
return self._velocity_obj.group_velocities[0, self._pp.band_indices, :]
def _set_gv_by_gv(self, i_gp, i_data):
@ -246,17 +369,42 @@ class ConductivityMixIn(HeatCapacityMixIn):
# Sum all vxv at k*
for j, vxv in enumerate(([0, 0], [1, 1], [2, 2], [1, 2], [0, 2], [0, 1])):
self._gv_sum2[i_data, :, j] = gv_by_gv_tensor[:, vxv[0], vxv[1]]
self._gv_by_gv[i_data, :, j] = gv_by_gv_tensor[:, vxv[0], vxv[1]]
def _get_gv_by_gv(self, i_gp, i_data):
multi = self._get_multiplicity_at_q(i_gp)
if self._is_kappa_star:
multi = get_multiplicity_at_q(
self._grid_points[i_gp], # type: ignore
self._pp,
self._point_operations,
)
else:
multi = 1
gv = self._gv[i_data]
gv_by_gv = np.zeros((len(gv), 3, 3), dtype="double")
for r in self._rotations_cartesian:
gvs_rot = np.dot(gv, r.T)
gv_by_gv += [np.outer(r_gv, r_gv) for r_gv in gvs_rot]
gv_by_gv /= multi
return gv_by_gv, self._get_kstar_order(i_gp, multi)
kstar_order = get_kstar_order(
self._grid_weights[i_gp], # type: ignore
multi,
self._point_operations,
verbose=self._log_level > 0,
)
return gv_by_gv, kstar_order
def _allocate_values(self):
super()._allocate_values()
num_band0 = len(self._pp.band_indices)
if self._is_reducible_collision_matrix:
num_grid_points = np.prod(self._pp.mesh_numbers)
else:
num_grid_points = len(self._grid_points)
self._gv_by_gv = np.zeros(
(num_grid_points, num_band0, 6), order="C", dtype="double"
)
class ConductivityBase(ABC):
@ -264,8 +412,6 @@ class ConductivityBase(ABC):
All Conductivity* classes have to inherit this base class.
self._gv has to be allocated in the inherited classes.
"""
_average_gv_over_kstar = False
@ -273,17 +419,16 @@ class ConductivityBase(ABC):
def __init__(
self,
interaction: Interaction,
grid_points=None,
temperatures: Optional[Union[List, np.ndarray]] = None,
sigmas: Optional[Union[List, np.ndarray]] = None,
sigma_cutoff: Optional[float] = None,
grid_points: np.ndarray | None = None,
temperatures: list | np.ndarray | None = None,
sigmas: list | np.ndarray | None = None,
sigma_cutoff: float | None = None,
is_isotope=False,
mass_variances: Optional[Union[List, np.ndarray]] = None,
boundary_mfp: Optional[float] = None,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
log_level=0,
mass_variances: list | np.ndarray | None = None,
boundary_mfp: float | None = None,
is_kappa_star: bool = True,
is_full_pp: bool = False,
log_level: int = 0,
):
"""Init method.
@ -319,10 +464,6 @@ class ConductivityBase(ABC):
iterating over specific grid points. With `is_kappa_star=True`
and `grid_points=None`, ir-grid points are used for the iteration.
Default is True.
gv_delta_q : float, optional, default is None, # for group velocity
With non-analytical correction, group velocity is calculated
by central finite difference method. This value gives the distance
in both directions in 1/Angstrom. The default value will be 1e-5.
is_full_pp : bool, optional, default is False
With True, full elements of phonon-phonon interaction strength
are computed. However with tetrahedron method, part of them are
@ -347,17 +488,14 @@ class ConductivityBase(ABC):
self._ir_grid_points,
self._grid_weights,
) = self._get_grid_info(grid_points)
self._grid_point_count: int = 0
self._num_sampling_grid_points: int = 0
self._grid_point_count = 0
self._sigmas: List
if sigmas is None:
self._sigmas = []
else:
self._sigmas = list(sigmas)
self._sigma_cutoff = sigma_cutoff
self._collision: Union[ImagSelfEnergy, CollisionMatrix]
self._temperatures: Optional[np.ndarray]
self._collision: ImagSelfEnergy | CollisionMatrix
if temperatures is None:
self._temperatures = None
else:
@ -370,7 +508,7 @@ class ConductivityBase(ABC):
self._eigenvectors,
self._phonon_done,
) = self._pp.get_phonons()
if (self._phonon_done == 0).any():
if not self._pp.phonon_all_done:
self._pp.run_phonon_solver()
self._is_isotope = is_isotope
@ -385,20 +523,13 @@ class ConductivityBase(ABC):
self._read_gamma_iso = False
# Allocated in self._allocate_values.
self._gv: np.ndarray
self._gamma: np.ndarray
self._gamma_iso: Optional[np.ndarray] = None
volume = self._pp.primitive.volume
self._conversion_factor = get_unit_to_WmK() / volume
self._conversion_factor = get_unit_to_WmK() / self._pp.primitive.volume
self._averaged_pp_interaction = None
# `self._velocity_obj` is the instance of an inherited class of
# `GroupVelocity`. `self._init_velocity()` is the method setup the instance,
# which must be implemented in the inherited class of `ConductivityBase`.
self._velocity_obj: GroupVelocity
self._init_velocity(gv_delta_q)
self._conductivity_components: ConductivityComponentsBase
def __iter__(self):
"""Calculate mode kappa at each grid point."""
@ -418,21 +549,29 @@ class ConductivityBase(ABC):
self._grid_point_count += 1
return self._grid_point_count - 1
@property
def mode_heat_capacities(self):
"""Return mode heat capacity at constant volume at grid points.
Grid points are those at mode kappa are calculated.
"""
return self._conductivity_components.mode_heat_capacities
@property
def group_velocities(self):
"""Return group velocities at grid points.
Grid points are those at mode kappa are calculated.
"""
return self._conductivity_components.group_velocities
@property
def mesh_numbers(self):
"""Return mesh numbers of GR-grid."""
return self._pp.mesh_numbers
def get_mesh_numbers(self):
"""Return mesh numbers of GR-grid."""
warnings.warn(
"Use attribute, Conductivity.mesh_numbers "
"instead of Conductivity.get_mesh_numbers().",
DeprecationWarning,
stacklevel=2,
)
return self.mesh_numbers
@property
def bz_grid(self):
"""Return GR-grid."""
@ -445,41 +584,18 @@ class ConductivityBase(ABC):
Grid points are those at mode kappa are calculated.
"""
assert self._frequencies is not None
return self._frequencies[self._grid_points]
def get_frequencies(self):
"""Return frequencies at grid points.
Grid points are those at mode kappa are calculated.
"""
warnings.warn(
"Use attribute, Conductivity.frequencies "
"instead of Conductivity.get_frequencies().",
DeprecationWarning,
stacklevel=2,
)
return self.frequencies
@property
def qpoints(self):
"""Return q-points where mode kappa are calculated."""
return np.array(
self._get_qpoint_from_gp_index(self._grid_points),
get_qpoints_from_bz_grid_points(self._grid_points, self._pp.bz_grid),
dtype="double",
order="C",
)
def get_qpoints(self):
"""Return q-points where mode kappa are calculated."""
warnings.warn(
"Use attribute, Conductivity.qpoints "
"instead of Conductivity.get_qpoints().",
DeprecationWarning,
stacklevel=2,
)
return self.qpoints
@property
def grid_points(self):
"""Return grid point indices where mode kappa are calculated.
@ -489,35 +605,11 @@ class ConductivityBase(ABC):
"""
return self._grid_points
def get_grid_points(self):
"""Return grid point indices where mode kappa are calculated.
Grid point indices are given in BZ-grid.
"""
warnings.warn(
"Use attribute, Conductivity.grid_points "
"instead of Conductivity.get_grid_points().",
DeprecationWarning,
stacklevel=2,
)
return self.grid_points
@property
def grid_weights(self):
"""Return grid point weights where mode kappa are calculated."""
return self._grid_weights
def get_grid_weights(self):
"""Return grid point weights where mode kappa are calculated."""
warnings.warn(
"Use attribute, Conductivity.grid_weights "
"instead of Conductivity.get_grid_weights().",
DeprecationWarning,
stacklevel=2,
)
return self.grid_weights
@property
def temperatures(self):
"""Setter and getter of temperatures."""
@ -528,26 +620,6 @@ class ConductivityBase(ABC):
self._temperatures = np.array(temperatures, dtype="double")
self._allocate_values()
def get_temperatures(self):
"""Return temperatures."""
warnings.warn(
"Use attribute, Conductivity.temperatures "
"instead of Conductivity.get_temperatures().",
DeprecationWarning,
stacklevel=2,
)
return self.temperatures
def set_temperatures(self, temperatures):
"""Set temperatures."""
warnings.warn(
"Use attribute, Conductivity.temperatures "
"instead of Conductivity.set_temperatures().",
DeprecationWarning,
stacklevel=2,
)
self.temperatures = temperatures
@property
def gamma(self):
"""Setter and getter of gamma."""
@ -558,24 +630,6 @@ class ConductivityBase(ABC):
self._gamma = gamma
self._read_gamma = True
def get_gamma(self):
"""Return gamma."""
warnings.warn(
"Use attribute, Conductivity.gamma instead of Conductivity.get_gamma().",
DeprecationWarning,
stacklevel=2,
)
return self.gamma
def set_gamma(self, gamma):
"""Set gamma."""
warnings.warn(
"Use attribute, Conductivity.gamma instead of Conductivity.set_gamma().",
DeprecationWarning,
stacklevel=2,
)
self.gamma = gamma
@property
def gamma_isotope(self):
"""Setter and getter of gamma from isotope."""
@ -586,100 +640,42 @@ class ConductivityBase(ABC):
self._gamma_iso = gamma_iso
self._read_gamma_iso = True
def get_gamma_isotope(self):
"""Return gamma from isotope."""
warnings.warn(
"Use attribute, Conductivity.gamma_isotope "
"instead of Conductivity.get_gamma_isotope().",
DeprecationWarning,
stacklevel=2,
)
return self.gamma_isotope
def set_gamma_isotope(self, gamma_iso):
"""Set gamma from isotope."""
warnings.warn(
"Use attribute, Conductivity.gamma_isotope "
"instead of Conductivity.set_gamma_isotope().",
DeprecationWarning,
stacklevel=2,
)
self.gamma_isotope = gamma_iso
@property
def sigmas(self):
"""Return sigmas."""
return self._sigmas
def get_sigmas(self):
"""Return sigmas."""
warnings.warn(
"Use attribute, Conductivity.sigmas instead of Conductivity.get_sigmas().",
DeprecationWarning,
stacklevel=2,
)
return self.sigmas
@property
def sigma_cutoff_width(self):
"""Return smearing width cutoff."""
return self._sigma_cutoff
def get_sigma_cutoff_width(self):
"""Return smearing width cutoff."""
warnings.warn(
"Use attribute, Conductivity.sigma_cutoff_width "
"instead of Conductivity.get_sigma_cutoff_width().",
DeprecationWarning,
stacklevel=2,
)
return self.sigma_cutoff_width
@property
def grid_point_count(self):
"""Return iterator count of self."""
return self._grid_point_count
def get_grid_point_count(self):
"""Return iterator count of self."""
warnings.warn(
"Use attribute, Conductivity.grid_point_count "
"instead of Conductivity.get_grid_point_count().",
DeprecationWarning,
stacklevel=2,
)
return self.grid_point_count
@property
def averaged_pp_interaction(self):
"""Return averaged pp strength."""
return self._averaged_pp_interaction
def get_averaged_pp_interaction(self):
"""Return averaged pp interaction strength."""
warnings.warn(
"Use attribute, Conductivity.averaged_pp_interaction "
"instead of Conductivity.get_averaged_pp_interaction().",
DeprecationWarning,
stacklevel=2,
)
return self.averaged_pp_interaction
@property
def boundary_mfp(self) -> float:
def boundary_mfp(self) -> Optional[float]:
"""Return boundary MFP."""
return self._boundary_mfp
def get_number_of_sampling_grid_points(self):
@property
def number_of_sampling_grid_points(self):
"""Return number of grid points.
This is calculated by the sum of numbers of arms of k-start in
`Conductivity._set_gv_by_gv`.
"""
return self._num_sampling_grid_points
return self._conductivity_components.number_of_sampling_grid_points
def _get_point_operations(self) -> Tuple[np.ndarray, np.ndarray]:
def _get_point_operations(self) -> tuple[np.ndarray, np.ndarray]:
"""Return reciprocal point group operations.
Returns
@ -705,7 +701,7 @@ class ConductivityBase(ABC):
return point_operations, rotations_cartesian
def _get_grid_info(self, grid_points) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
def _get_grid_info(self, grid_points) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Return grid point information in BZGrid.
Returns
@ -738,47 +734,22 @@ class ConductivityBase(ABC):
@abstractmethod
def _run_at_grid_point(self):
"""Run at conductivity calculation at specified grid point.
Should be implemented in Conductivity* class.
"""
"""Run at conductivity calculation at specified grid point."""
raise NotImplementedError()
@abstractmethod
def _allocate_values(self):
"""Allocate necessary data arrays.
Should be implemented in Conductivity* class.
"""
"""Allocate necessary data arrays."""
raise NotImplementedError()
@abstractmethod
def _set_velocities(self, i_gp, i_data):
"""Set velocities at grid point and at data location.
Should be implemented in Conductivity*MixIn.
"""
raise NotImplementedError()
@abstractmethod
def _init_velocity(self, gv_delta_q):
"""Initialize velocity class instance.
Should be implemented in Conductivity*MixIn.
"""
"""Set velocities at grid point and at data location."""
raise NotImplementedError()
@abstractmethod
def _set_cv(self, i_gp, i_data):
"""Set heat capacity at grid point and at data location.
Should be implemented in Conductivity*MixIn.
"""
"""Set heat capacity at grid point and at data location."""
raise NotImplementedError()
def _get_ir_grid_points(self, grid_points):
@ -802,60 +773,6 @@ class ConductivityBase(ABC):
return ir_grid_points, grid_weights
def _get_qpoint_from_gp_index(self, i_gps):
"""Return q-point(s) in reduced coordinates of grid point(s).
Parameters
----------
i_gps : int or ndarray
BZ-grid index (int) or indices (ndarray).
"""
return np.dot(self._pp.bz_grid.addresses[i_gps], self._pp.bz_grid.QDinv.T)
def _get_multiplicity_at_q(self, i_gp):
"""Return multiplicity (order of site-symmetry) of q-point."""
if self._is_kappa_star:
q = self._get_qpoint_from_gp_index(self._grid_points[i_gp])
reclat = np.linalg.inv(self._pp.primitive.cell)
multi = 0
for q_rot in [np.dot(r, q) for r in self._point_operations]:
diff = q - q_rot
diff -= np.rint(diff)
dist = np.linalg.norm(np.dot(reclat, diff))
if dist < self._pp.primitive_symmetry.tolerance:
multi += 1
else:
multi = 1
return multi
def _get_kstar_order(self, i_gp, multi):
"""Return order (number of arms) of kstar.
multi : int
Multiplicity of q-point of `i_gp`, which can be obtained by
`self._get_multiplicity_at_q(i_gp)`.
"""
order_kstar = len(self._point_operations) // multi
if order_kstar != self._grid_weights[i_gp]:
if self._log_level:
text = (
"Number of elements in k* is unequal "
"to number of equivalent grid-points. "
"This means that the mesh sampling grids break "
"symmetry. Please check carefully "
"the convergence over grid point densities."
)
msg = textwrap.fill(
text, initial_indent=" ", subsequent_indent=" ", width=70
)
print("*" * 30 + "Warning" + "*" * 30)
print(msg)
print("*" * 67)
return order_kstar
def _get_gamma_isotope_at_sigmas(self, i):
gamma_iso = []
for sigma in self._sigmas:
@ -906,7 +823,7 @@ class ConductivityBase(ABC):
num_band = len(self._pp.primitive) * 3
g_boundary = np.zeros(num_band, dtype="double")
try:
gv = self._gv
gv = self._conductivity_components.group_velocities
except AttributeError:
print("(_get_boundary_scattering) _gv has to be implemented.")
return g_boundary
@ -927,7 +844,7 @@ class ConductivityBase(ABC):
"======================= Grid point %d (%d/%d) "
"=======================" % (bzgp, i_gp + 1, len(self._grid_points))
)
qpoint = self._get_qpoint_from_gp_index(bzgp)
qpoint = get_qpoints_from_bz_grid_points(bzgp, self._pp.bz_grid)
print("q-point: (%5.2f %5.2f %5.2f)" % tuple(qpoint))
if self._boundary_mfp is not None:
if self._boundary_mfp > 1000:

107
phono3py/conductivity/direct_solution.py
View File

@ -34,11 +34,13 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import warnings
from __future__ import annotations
from typing import Optional, Union
import numpy as np
from phono3py.conductivity.base import ConductivityMixIn
from phono3py.conductivity.base import ConductivityComponents
from phono3py.conductivity.direct_solution_base import (
ConductivityLBTEBase,
diagonalize_collision_matrix,
@ -46,38 +48,38 @@ from phono3py.conductivity.direct_solution_base import (
from phono3py.phonon3.interaction import Interaction
class ConductivityLBTE(ConductivityMixIn, ConductivityLBTEBase):
class ConductivityLBTE(ConductivityLBTEBase):
"""Lattice thermal conductivity calculation by direct solution."""
def __init__(
self,
interaction: Interaction,
grid_points=None,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
is_isotope=False,
mass_variances=None,
boundary_mfp=None,
solve_collective_phonon=False,
is_reducible_collision_matrix=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
read_pp=False,
pp_filename=None,
pinv_cutoff=1.0e-8,
pinv_solver=0,
pinv_method=0,
log_level=0,
lang="C",
grid_points: Optional[np.ndarray] = None,
temperatures: Optional[Union[list, np.ndarray]] = None,
sigmas: Optional[Union[list, np.ndarray]] = None,
sigma_cutoff: Optional[float] = None,
is_isotope: bool = False,
mass_variances: Optional[Union[list, np.ndarray]] = None,
boundary_mfp: Optional[float] = None, # in micrometer
solve_collective_phonon: bool = False,
is_reducible_collision_matrix: bool = False,
is_kappa_star: bool = True,
gv_delta_q: Optional[float] = None,
is_full_pp: bool = False,
read_pp: bool = False,
pp_filename: Optional[float] = None,
pinv_cutoff: float = 1.0e-8,
pinv_solver: int = 0,
pinv_method: int = 0,
log_level: int = 0,
lang: str = "C",
):
"""Init method."""
self._kappa = None
self._kappa_RTA = None
self._mode_kappa = None
self._mode_kappa_RTA = None
self._gv_sum2 = None
self._gv_by_gv = None
super().__init__(
interaction,
@ -91,7 +93,6 @@ class ConductivityLBTE(ConductivityMixIn, ConductivityLBTEBase):
solve_collective_phonon=solve_collective_phonon,
is_reducible_collision_matrix=is_reducible_collision_matrix,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_full_pp=is_full_pp,
read_pp=read_pp,
pp_filename=pp_filename,
@ -102,35 +103,47 @@ class ConductivityLBTE(ConductivityMixIn, ConductivityLBTEBase):
lang=lang,
)
self._conductivity_components = ConductivityComponents(
self._pp,
self._grid_points,
self._grid_weights,
self._point_operations,
self._rotations_cartesian,
temperatures=self._temperatures,
average_gv_over_kstar=self._average_gv_over_kstar,
is_kappa_star=self._is_kappa_star,
gv_delta_q=gv_delta_q,
is_reducible_collision_matrix=self._is_reducible_collision_matrix,
log_level=self._log_level,
)
@property
def kappa(self):
"""Return kappa."""
return self._kappa
@property
def mode_kappa(self):
"""Return mode_kappa."""
return self._mode_kappa
@property
def kappa_RTA(self):
"""Return RTA lattice thermal conductivity."""
return self._kappa_RTA
def get_kappa_RTA(self):
"""Return RTA lattice thermal conductivity."""
warnings.warn(
"Use attribute, Conductivity_LBTE.kappa_RTA "
"instead of Conductivity_LBTE.get_kappa_RTA().",
DeprecationWarning,
stacklevel=2,
)
return self.kappa_RTA
@property
def mode_kappa_RTA(self):
"""Return RTA mode lattice thermal conductivities."""
return self._mode_kappa_RTA
def get_mode_kappa_RTA(self):
"""Return RTA mode lattice thermal conductivities."""
warnings.warn(
"Use attribute, Conductivity_LBTE.mode_kappa_RTA "
"instead of Conductivity_LBTE.get_mode_kappa_RTA().",
DeprecationWarning,
stacklevel=2,
)
return self.mode_kappa_RTA
def _set_cv(self, i_gp, i_data):
"""Set cv for conductivity components."""
self._conductivity_components.set_heat_capacities(i_gp, i_data)
def _set_velocities(self, i_gp, i_data):
"""Set velocities for conductivity components."""
self._conductivity_components.set_velocities(i_gp, i_data)
def _allocate_local_values(self, num_grid_points):
"""Allocate grid point local arrays.
@ -140,6 +153,12 @@ class ConductivityLBTE(ConductivityMixIn, ConductivityLBTEBase):
grid points to be iterated over.
"""
if self._temperatures is None:
raise RuntimeError(
"Temperatures have not been set yet. "
"Set temperatures before this method."
)
num_band0 = len(self._pp.band_indices)
num_temp = len(self._temperatures)
super()._allocate_local_values(num_grid_points)
@ -150,7 +169,7 @@ class ConductivityLBTE(ConductivityMixIn, ConductivityLBTEBase):
self._kappa_RTA = np.zeros(
(len(self._sigmas), num_temp, 6), dtype="double", order="C"
)
self._gv_sum2 = np.zeros(
self._gv_by_gv = np.zeros(
(num_grid_points, num_band0, 6), dtype="double", order="C"
)
self._mode_kappa = np.zeros(

122
phono3py/conductivity/direct_solution_base.py
View File

@ -34,11 +34,12 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import sys
import time
import warnings
from abc import abstractmethod
from typing import Optional
from typing import Optional, Union
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
@ -62,25 +63,24 @@ class ConductivityLBTEBase(ConductivityBase):
def __init__(
self,
interaction: Interaction,
grid_points=None,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
is_isotope=False,
mass_variances=None,
boundary_mfp=None, # in micrometer
solve_collective_phonon=False,
is_reducible_collision_matrix=False,
is_kappa_star=True,
gv_delta_q=None, # finite difference for group velocity
is_full_pp=False,
read_pp=False,
pp_filename=None,
pinv_cutoff=1.0e-8,
pinv_solver=0,
pinv_method=0,
log_level=0,
lang="C",
grid_points: Optional[np.ndarray] = None,
temperatures: Optional[Union[list, np.ndarray]] = None,
sigmas: Optional[Union[list, np.ndarray]] = None,
sigma_cutoff: Optional[float] = None,
is_isotope: bool = False,
mass_variances: Optional[Union[list, np.ndarray]] = None,
boundary_mfp: Optional[float] = None, # in micrometer
solve_collective_phonon: bool = False,
is_reducible_collision_matrix: bool = False,
is_kappa_star: bool = True,
is_full_pp: bool = False,
read_pp: bool = False,
pp_filename: Optional[float] = None,
pinv_cutoff: float = 1.0e-8,
pinv_solver: int = 0,
pinv_method: int = 0,
log_level: int = 0,
lang: str = "C",
):
"""Init method."""
super().__init__(
@ -97,8 +97,6 @@ class ConductivityLBTEBase(ConductivityBase):
log_level=log_level,
)
self._init_velocity(gv_delta_q)
self._lang = lang
self._collision_eigenvalues = None
self._is_reducible_collision_matrix = is_reducible_collision_matrix
@ -150,41 +148,11 @@ class ConductivityLBTEBase(ConductivityBase):
def collision_matrix(self, collision_matrix):
self._collision_matrix = collision_matrix
def get_collision_matrix(self):
"""Return collision matrix."""
warnings.warn(
"Use attribute, Conductivity_LBTE.collision_matrix "
"instead of Conductivity_LBTE.get_collision_matrix().",
DeprecationWarning,
stacklevel=2,
)
return self.collision_matrix
def set_collision_matrix(self, collision_matrix):
"""Set collision matrix."""
warnings.warn(
"Use attribute, Conductivity_LBTE.collision_matrix "
"instead of Conductivity_LBTE.set_collision_matrix().",
DeprecationWarning,
stacklevel=2,
)
self.collision_matrix = collision_matrix
@property
def collision_eigenvalues(self):
"""Return eigenvalues of collision matrix."""
return self._collision_eigenvalues
def get_collision_eigenvalues(self):
"""Return eigenvalues of collision matrix."""
warnings.warn(
"Use attribute, Conductivity_LBTE.collision_eigenvalues "
"instead of Conductivity_LBTE.get_collision_eigenvalues().",
DeprecationWarning,
stacklevel=2,
)
return self.collision_eigenvalues
def get_frequencies_all(self):
"""Return phonon frequencies on GR-grid."""
return self._frequencies[self._pp.bz_grid.grg2bzg]
@ -239,7 +207,7 @@ class ConductivityLBTEBase(ConductivityBase):
kappa and mode_kappa are overwritten.
"""
N = self._num_sampling_grid_points
N = self.number_of_sampling_grid_points
if self._solve_collective_phonon:
self._set_mode_kappa_Chaput(mode_kappa, i_sigma, i_temp, weights)
else:
@ -269,7 +237,7 @@ class ConductivityLBTEBase(ConductivityBase):
kappa and mode_kappa are overwritten.
"""
N = self._num_sampling_grid_points
N = self.number_of_sampling_grid_points
X = self._get_X(i_temp, weights)
num_mesh_points = np.prod(self._pp.mesh_numbers)
Y = self._get_Y(i_sigma, i_temp, weights, X)
@ -314,10 +282,6 @@ class ConductivityLBTEBase(ConductivityBase):
"""Allocate grid point local arrays."""
num_band0 = len(self._pp.band_indices)
num_temp = len(self._temperatures)
self._gv = np.zeros((num_grid_points, num_band0, 3), dtype="double", order="C")
self._cv = np.zeros(
(num_temp, num_grid_points, num_band0), dtype="double", order="C"
)
self._gamma = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band0),
dtype="double",
@ -625,12 +589,20 @@ class ConductivityLBTEBase(ConductivityBase):
sys.stdout.flush()
def _expand_local_values(self, ir_gr_grid_points, rot_grid_points):
"""Fill elements of local properties at grid points."""
"""Fill elements of local properties at grid points.
Note
----
Internal state of self._conductivity_components is updated.
"""
cv = self._conductivity_components.mode_heat_capacities
gv = self._conductivity_components.group_velocities
for ir_gp in ir_gr_grid_points:
cv_irgp = self._cv[:, ir_gp, :].copy()
self._cv[:, ir_gp, :] = 0
gv_irgp = self._gv[ir_gp].copy()
self._gv[ir_gp] = 0
cv_irgp = cv[:, ir_gp, :].copy()
cv[:, ir_gp, :] = 0
gv_irgp = gv[ir_gp].copy()
gv[ir_gp] = 0
gamma_irgp = self._gamma[:, :, ir_gp, :].copy()
self._gamma[:, :, ir_gp, :] = 0
multi = (rot_grid_points[:, ir_gp] == ir_gp).sum()
@ -642,8 +614,8 @@ class ConductivityLBTEBase(ConductivityBase):
self._gamma[:, :, gp_r, :] += gamma_irgp / multi
if self._is_isotope:
self._gamma_iso[:, gp_r, :] += gamma_iso_irgp / multi
self._cv[:, gp_r, :] += cv_irgp / multi
self._gv[gp_r] += np.dot(gv_irgp, r.T) / multi
cv[:, gp_r, :] += cv_irgp / multi
gv[gp_r] += np.dot(gv_irgp, r.T) / multi
def _get_weights(self):
"""Return weights used for collision matrix and |X> and |f>.
@ -782,7 +754,7 @@ class ConductivityLBTEBase(ConductivityBase):
def _get_X(self, i_temp, weights):
"""Calculate X in Chaput's paper."""
num_band = len(self._pp.primitive) * 3
X = self._gv.copy()
X = self._conductivity_components.group_velocities.copy()
if self._is_reducible_collision_matrix:
freqs = self._frequencies[self._pp.bz_grid.grg2bzg]
else:
@ -963,7 +935,7 @@ class ConductivityLBTEBase(ConductivityBase):
This RTA is supposed to be the same as conductivity_RTA.
"""
N = self._num_sampling_grid_points
N = self.number_of_sampling_grid_points
num_band = len(self._pp.primitive) * 3
X = self._get_X(i_temp, weights)
Y = np.zeros_like(X)
@ -1013,7 +985,7 @@ class ConductivityLBTEBase(ConductivityBase):
`kappa` and `mode_kappa` are overwritten.
"""
N = self._num_sampling_grid_points
N = self.number_of_sampling_grid_points
num_band = len(self._pp.primitive) * 3
X = self._get_X(i_temp, weights)
Y = np.zeros_like(X)
@ -1156,7 +1128,7 @@ class ConductivityLBTEBase(ConductivityBase):
# shape = (num_grid_points, num_band, 3),
for i, f_gp in enumerate(self._f_vectors):
for j, f in enumerate(f_gp):
cv = self._cv[i_temp, i, j]
cv = self._conductivity_components.mode_heat_capacities[i_temp, i, j]
if cv < 1e-10:
continue
self._mfp[i_sigma, i_temp, i, j] = (
@ -1167,19 +1139,21 @@ class ConductivityLBTEBase(ConductivityBase):
gp = self._grid_points[i]
frequencies = self._frequencies[gp]
if self._is_reducible_collision_matrix:
gv = self._gv[self._pp.bz_grid.bzg2grg[gp]]
gv = self._conductivity_components.group_velocities[
self._pp.bz_grid.bzg2grg[gp]
]
else:
gv = self._gv[i]
gv = self._conductivity_components.group_velocities[i]
if self._is_full_pp:
ave_pp = self._averaged_pp_interaction[i]
text = "Frequency group velocity (x, y, z) |gv| Pqj"
else:
text = "Frequency group velocity (x, y, z) |gv|"
if self._velocity_obj.q_length is None:
if self._conductivity_components.gv_delta_q is None:
pass
else:
text += " (dq=%3.1e)" % self._velocity_obj.q_length
text += " (dq=%3.1e)" % self._conductivity_components.gv_delta_q
print(text)
if self._is_full_pp:
for f, v, pp in zip(frequencies, gv, ave_pp):

4
phono3py/conductivity/init_rta.py
View File

@ -182,7 +182,7 @@ class ShowCalcProgress:
kappa = br.kappa
num_ignored_phonon_modes = br.number_of_ignored_phonon_modes
num_band = br.frequencies.shape[1]
num_phonon_modes = br.get_number_of_sampling_grid_points() * num_band
num_phonon_modes = br.number_of_sampling_grid_points * num_band
for i, sigma in enumerate(sigmas):
text = "----------- Thermal conductivity (W/m-k) "
if sigma:
@ -223,7 +223,7 @@ class ShowCalcProgress:
kappa_C = br.kappa_C
num_ignored_phonon_modes = br.number_of_ignored_phonon_modes
num_band = br.frequencies.shape[1]
num_phonon_modes = br.get_number_of_sampling_grid_points() * num_band
num_phonon_modes = br.number_of_sampling_grid_points * num_band
for i, sigma in enumerate(sigmas):
text = "----------- Thermal conductivity (W/m-k) "
if sigma:

26
phono3py/conductivity/kubo_base.py
View File

@ -37,6 +37,8 @@
import numpy as np
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.base import get_kstar_order, get_multiplicity_at_q
from phono3py.phonon.grid import get_qpoints_from_bz_grid_points
from phono3py.phonon.group_velocity_matrix import GroupVelocityMatrix
from phono3py.phonon.heat_capacity_matrix import mode_cv_matrix
@ -102,7 +104,9 @@ class ConductivityKuboMixIn:
"""
irgp = self._grid_points[i_gp]
self._velocity_obj.run([self._get_qpoint_from_gp_index(irgp)])
self._velocity_obj.run(
[get_qpoints_from_bz_grid_points(irgp, self._pp.bz_grid)]
)
gvm = np.zeros(self._gv_mat.shape[1:], dtype=self._complex_dtype, order="C")
gv = np.zeros(self._gv.shape[1:], dtype="double", order="C")
for i in range(3):
@ -136,8 +140,15 @@ class ConductivityKuboMixIn:
Number of kstar arms.
"""
multi = self._get_multiplicity_at_q(i_gp)
q = self._get_qpoint_from_gp_index(self._grid_points[i_gp])
if self._is_kappa_star:
multi = get_multiplicity_at_q(
self._grid_points[i_gp],
self._pp,
self._point_operations,
)
else:
multi = 1
q = get_qpoints_from_bz_grid_points(self._grid_points[i_gp], self._pp.bz_grid)
qpoints = [np.dot(r, q) for r in self._point_operations]
self._velocity_obj.run(qpoints)
@ -151,5 +162,10 @@ class ConductivityKuboMixIn:
gvm_by_gvm = np.multiply(gvm[a], gvm[b].T)
gvm_sum2[:, :, i_pair] += gvm_by_gvm[self._pp.band_indices, :]
gvm_sum2 /= multi
return gvm_sum2, self._get_kstar_order(i_gp, multi)
kstar_order = get_kstar_order(
self._grid_weights[i_gp],
multi,
self._point_operations,
verbose=self._log_level > 0,
)
return gvm_sum2, kstar_order

122
phono3py/conductivity/rta.py
View File

@ -34,43 +34,45 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
from typing import Optional, Union
import numpy as np
from phono3py.conductivity.base import ConductivityMixIn
from phono3py.conductivity.base import ConductivityComponents
from phono3py.conductivity.rta_base import ConductivityRTABase
from phono3py.phonon3.interaction import Interaction
class ConductivityRTA(ConductivityMixIn, ConductivityRTABase):
class ConductivityRTA(ConductivityRTABase):
"""Lattice thermal conductivity calculation with RTA."""
def __init__(
self,
interaction: Interaction,
grid_points=None,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
is_isotope=False,
mass_variances=None,
boundary_mfp=None, # in micrometer
use_ave_pp=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
read_pp=False,
store_pp=False,
pp_filename=None,
is_N_U=False,
is_gamma_detail=False,
is_frequency_shift_by_bubble=False,
log_level=0,
grid_points: Optional[np.ndarray] = None,
temperatures: Optional[Union[list, np.ndarray]] = None,
sigmas: Optional[Union[list, np.ndarray]] = None,
sigma_cutoff: Optional[float] = None,
is_isotope: bool = False,
mass_variances: Optional[Union[list, np.ndarray]] = None,
boundary_mfp: Optional[float] = None, # in micrometer
use_ave_pp: bool = False,
is_kappa_star: bool = True,
gv_delta_q: Optional[float] = None,
is_full_pp: bool = False,
read_pp: bool = False,
store_pp: bool = False,
pp_filename: Optional[float] = None,
is_N_U: bool = False,
is_gamma_detail: bool = False,
is_frequency_shift_by_bubble: bool = False,
log_level: int = 0,
):
"""Init method."""
self._cv = None
self._kappa = None
self._mode_kappa = None
self._gv_sum2 = None
super().__init__(
interaction,
@ -83,7 +85,6 @@ class ConductivityRTA(ConductivityMixIn, ConductivityRTABase):
boundary_mfp=boundary_mfp,
use_ave_pp=use_ave_pp,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_full_pp=is_full_pp,
read_pp=read_pp,
store_pp=store_pp,
@ -94,13 +95,62 @@ class ConductivityRTA(ConductivityMixIn, ConductivityRTABase):
log_level=log_level,
)
self._conductivity_components = ConductivityComponents(
self._pp,
self._grid_points,
self._grid_weights,
self._point_operations,
self._rotations_cartesian,
temperatures=self._temperatures,
average_gv_over_kstar=self._average_gv_over_kstar,
is_kappa_star=self._is_kappa_star,
gv_delta_q=gv_delta_q,
log_level=self._log_level,
)
@property
def kappa(self):
"""Return kappa."""
return self._kappa
@property
def mode_kappa(self):
"""Return mode_kappa."""
return self._mode_kappa
@property
def gv_by_gv(self):
"""Return gv_by_gv at grid points where mode kappa are calculated."""
return self._conductivity_components.gv_by_gv
def _set_cv(self, i_gp, i_data):
"""Set cv for conductivity components."""
self._conductivity_components.set_heat_capacities(i_gp, i_data)
def _set_velocities(self, i_gp, i_data):
"""Set velocities for conductivity components."""
self._conductivity_components.set_velocities(i_gp, i_data)
def set_kappa_at_sigmas(self):
"""Calculate kappa from ph-ph interaction results."""
if not self._pp.phonon_all_done:
raise RuntimeError(
"Phonon calculation has not been done yet. "
"Run phono3py.run_phonon_solver() before this method."
)
if self._temperatures is None:
raise RuntimeError(
"Temperatures have not been set yet. "
"Set temperatures before this method."
)
num_band = len(self._pp.primitive) * 3
mode_heat_capacities = self._conductivity_components.mode_heat_capacities
gv_by_gv = self._conductivity_components.gv_by_gv
for i, _ in enumerate(self._grid_points):
cv = self._cv[:, i, :]
cv = mode_heat_capacities[:, i, :]
gp = self._grid_points[i]
frequencies = self._frequencies[gp]
frequencies = self._frequencies[gp] # type: ignore
# Kappa
for j in range(len(self._sigmas)):
@ -108,13 +158,13 @@ class ConductivityRTA(ConductivityMixIn, ConductivityRTABase):
g_sum = self._get_main_diagonal(i, j, k)
for ll in range(num_band):
if frequencies[ll] < self._pp.cutoff_frequency:
self._num_ignored_phonon_modes[j, k] += 1
self._num_ignored_phonon_modes[j, k] += 1 # type: ignore
continue
old_settings = np.seterr(all="raise")
try:
self._mode_kappa[j, k, i, ll] = (
self._gv_sum2[i, ll]
self._mode_kappa[j, k, i, ll] = ( # type: ignore
gv_by_gv[i, ll]
* cv[k, ll]
/ (g_sum[ll] * 2)
* self._conversion_factor
@ -135,24 +185,22 @@ class ConductivityRTA(ConductivityMixIn, ConductivityRTABase):
print("=" * 61)
np.seterr(**old_settings)
N = self._num_sampling_grid_points
N = self.number_of_sampling_grid_points
self._kappa = self._mode_kappa.sum(axis=2).sum(axis=2) / N
def _allocate_values(self):
if self._temperatures is None:
raise RuntimeError(
"Temperatures have not been set yet. "
"Set temperatures before this method."
)
super()._allocate_values()
num_band0 = len(self._pp.band_indices)
num_band = len(self._pp.primitive) * 3
num_grid_points = len(self._grid_points)
num_temp = len(self._temperatures)
self._cv = np.zeros(
(num_temp, num_grid_points, num_band0), order="C", dtype="double"
)
self._gv_sum2 = np.zeros(
(num_grid_points, num_band0, 6), order="C", dtype="double"
)
# kappa* and mode_kappa* are accessed when all bands exist, i.e.,
# num_band0==num_band.
self._kappa = np.zeros(

65
phono3py/conductivity/rta_base.py
View File

@ -34,6 +34,8 @@
# 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 abc import abstractmethod
@ -43,7 +45,10 @@ from phonopy.physical_units import get_physical_units
from phono3py.conductivity.base import ConductivityBase
from phono3py.file_IO import read_pp_from_hdf5
from phono3py.other.tetrahedron_method import get_tetrahedra_relative_grid_address
from phono3py.phonon.grid import get_grid_points_by_rotations
from phono3py.phonon.grid import (
get_grid_points_by_rotations,
get_qpoints_from_bz_grid_points,
)
from phono3py.phonon3.imag_self_energy import ImagSelfEnergy, average_by_degeneracy
from phono3py.phonon3.interaction import Interaction
@ -58,24 +63,23 @@ class ConductivityRTABase(ConductivityBase):
def __init__(
self,
interaction: Interaction,
grid_points=None,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
is_isotope=False,
mass_variances=None,
boundary_mfp=None, # in micrometer
use_ave_pp=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
read_pp=False,
store_pp=False,
pp_filename=None,
is_N_U=False,
is_gamma_detail=False,
is_frequency_shift_by_bubble=False,
log_level=0,
grid_points: np.ndarray | None = None,
temperatures: list | np.ndarray | None = None,
sigmas: list | np.ndarray | None = None,
sigma_cutoff: float | None = None,
is_isotope: bool = False,
mass_variances: list | np.ndarray | None = None,
boundary_mfp: float | None = None, # in micrometer
use_ave_pp: bool = False,
is_kappa_star: bool = True,
is_full_pp: bool = False,
read_pp: bool = False,
store_pp: bool = False,
pp_filename: float | None = None,
is_N_U: bool = False,
is_gamma_detail: bool = False,
is_frequency_shift_by_bubble: bool = False,
log_level: int = 0,
):
"""Init method."""
self._is_N_U = is_N_U
@ -99,7 +103,6 @@ class ConductivityRTABase(ConductivityBase):
mass_variances=mass_variances,
boundary_mfp=boundary_mfp,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_full_pp=is_full_pp,
log_level=log_level,
)
@ -153,6 +156,12 @@ class ConductivityRTABase(ConductivityBase):
raise NotImplementedError()
def _allocate_values(self):
if self._temperatures is None:
raise RuntimeError(
"Temperatures have not been set yet. "
"Set temperatures before this method."
)
num_band0 = len(self._pp.band_indices)
num_grid_points = len(self._grid_points)
num_temp = len(self._temperatures)
@ -166,7 +175,6 @@ class ConductivityRTABase(ConductivityBase):
self._gamma_N = np.zeros_like(self._gamma)
self._gamma_U = np.zeros_like(self._gamma)
self._gv = np.zeros((num_grid_points, num_band0, 3), order="C", dtype="double")
if self._is_isotope:
self._gamma_iso = np.zeros(
(len(self._sigmas), num_grid_points, num_band0),
@ -320,7 +328,6 @@ class ConductivityRTABase(ConductivityBase):
s2p,
masses,
) = self._pp.get_primitive_and_supercell_correspondence()
fc3 = self._pp.fc3
triplets_at_q, weights_at_q, _, _ = self._pp.get_triplets_at_q()
if None in self._sigmas:
@ -376,7 +383,8 @@ class ConductivityRTABase(ConductivityBase):
self._pp.bz_grid.store_dense_gp_map * 1 + 1,
self._pp.bz_grid.D_diag,
self._pp.bz_grid.Q,
fc3,
self._pp.fc3,
self._pp.fc3_nonzero_indices,
svecs,
multi,
masses,
@ -407,7 +415,8 @@ class ConductivityRTABase(ConductivityBase):
self._pp.bz_grid.addresses,
self._pp.bz_grid.D_diag,
self._pp.bz_grid.Q,
fc3,
self._pp.fc3,
self._pp.fc3_nonzero_indices,
svecs,
multi,
masses,
@ -449,10 +458,10 @@ class ConductivityRTABase(ConductivityBase):
)
def _show_log(self, i_gp):
q = self._get_qpoint_from_gp_index(i_gp)
q = get_qpoints_from_bz_grid_points(i_gp, self._pp.bz_grid)
gp = self._grid_points[i_gp]
frequencies = self._frequencies[gp][self._pp.band_indices]
gv = self._gv[i_gp]
gv = self._conductivity_components.group_velocities[i_gp]
if self._averaged_pp_interaction is not None:
ave_pp = self._averaged_pp_interaction[i_gp]
@ -512,8 +521,8 @@ class ConductivityRTABase(ConductivityBase):
text = "Frequency group velocity (x, y, z) |gv| Pqj"
else:
text = "Frequency group velocity (x, y, z) |gv|"
if self._velocity_obj.q_length is None:
if self._conductivity_components.gv_delta_q is None:
pass
else:
text += " (dq=%3.1e)" % self._velocity_obj.q_length
text += " (dq=%3.1e)" % self._conductivity_components.gv_delta_q
print(text)

147
phono3py/conductivity/wigner_base.py
View File

@ -35,47 +35,83 @@
# POSSIBILITY OF SUCH DAMAGE.
import textwrap
from typing import Optional
import numpy as np
from numpy.typing import NDArray
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.base import HeatCapacityMixIn
from phono3py.phonon.grid import get_grid_points_by_rotations
from phono3py.conductivity.base import ConductivityComponentsBase, get_heat_capacities
from phono3py.phonon.grid import (
get_grid_points_by_rotations,
get_qpoints_from_bz_grid_points,
)
from phono3py.phonon.velocity_operator import VelocityOperator
from phono3py.phonon3.interaction import Interaction
class ConductivityWignerMixIn(HeatCapacityMixIn):
"""Thermal conductivity mix-in for velocity operator.
def get_conversion_factor_WTE(volume):
"""Return conversion factor of thermal conductivity."""
return (
(get_physical_units().THz * get_physical_units().Angstrom)
** 2 # --> group velocity
* get_physical_units().EV # --> specific heat is in eV/
* get_physical_units().Hbar # --> transform lorentzian_div_hbar from eV^-1 to s
/ (volume * get_physical_units().Angstrom ** 3)
) # --> unit cell volume
This mix-in is included in `ConductivityWignerRTA` and `ConductivityWignerLBTE`.
class ConductivityWignerComponents(ConductivityComponentsBase):
"""Thermal conductivity components for velocity operator.
Used by `ConductivityWignerRTA` and `ConductivityWignerLBTE`.
"""
@property
def kappa_TOT_RTA(self):
"""Return kappa."""
return self._kappa_TOT_RTA
def __init__(
self,
pp: Interaction,
grid_points: NDArray[np.int64],
grid_weights: NDArray[np.int64],
point_operations: NDArray[np.int64],
rotations_cartesian: NDArray[np.int64],
temperatures: Optional[NDArray[np.float64]] = None,
is_kappa_star: bool = True,
gv_delta_q: Optional[float] = None,
is_reducible_collision_matrix: bool = False,
log_level: int = 0,
):
"""Init method."""
super().__init__(
pp,
grid_points=grid_points,
grid_weights=grid_weights,
point_operations=point_operations,
rotations_cartesian=rotations_cartesian,
temperatures=temperatures,
is_kappa_star=is_kappa_star,
is_reducible_collision_matrix=is_reducible_collision_matrix,
log_level=log_level,
)
@property
def kappa_P_RTA(self):
"""Return kappa."""
return self._kappa_P_RTA
self._gv_operator: np.ndarray
self._gv_operator_sum2: np.ndarray
@property
def kappa_C(self):
"""Return kappa."""
return self._kappa_C
if self._pp.dynamical_matrix is None:
raise RuntimeError("Interaction.init_dynamical_matrix() has to be called.")
self._velocity_obj = VelocityOperator(
self._pp.dynamical_matrix,
q_length=gv_delta_q,
symmetry=self._pp.primitive_symmetry,
frequency_factor_to_THz=self._pp.frequency_factor_to_THz,
)
@property
def mode_kappa_P_RTA(self):
"""Return mode_kappa."""
return self._mode_kappa_P_RTA
self._num_sampling_grid_points = 0
self._complex_dtype = "c%d" % (np.dtype("double").itemsize * 2)
@property
def mode_kappa_C(self):
"""Return mode_kappa."""
return self._mode_kappa_C
if self._temperatures is not None:
self._allocate_values()
@property
def velocity_operator(self):
@ -91,27 +127,34 @@ class ConductivityWignerMixIn(HeatCapacityMixIn):
"""Return gv_by_gv operator at grid points where mode kappa are calculated."""
return self._gv_operator_sum2
def _init_velocity(self, gv_delta_q):
self._velocity_obj = VelocityOperator(
self._pp.dynamical_matrix,
q_length=gv_delta_q,
symmetry=self._pp.primitive_symmetry,
frequency_factor_to_THz=self._pp.frequency_factor_to_THz,
)
def _set_velocities(self, i_gp, i_data):
def set_velocities(self, i_gp, i_data):
"""Set velocities at a grid point."""
self._set_gv_operator(i_gp, i_data)
self._set_gv_by_gv_operator(i_gp, i_data)
def set_heat_capacities(self, i_gp: int, i_data: int):
"""Set heat capacity at grid point and at data location."""
if self._temperatures is None:
raise RuntimeError(
"Temperatures have not been set yet. "
"Set temperatures before this method."
)
cv = get_heat_capacities(self._grid_points[i_gp], self._pp, self._temperatures)
self._cv[:, i_data, :] = cv
def _set_gv_operator(self, i_irgp, i_data):
"""Set velocity operator."""
irgp = self._grid_points[i_irgp]
self._velocity_obj.run([self._get_qpoint_from_gp_index(irgp)])
frequencies, _, _ = self._pp.get_phonons()
self._velocity_obj.run(
[get_qpoints_from_bz_grid_points(irgp, self._pp.bz_grid)]
)
gv_operator = self._velocity_obj.velocity_operators[0, :, :, :]
self._gv_operator[i_data] = gv_operator[self._pp.band_indices, :, :]
#
gv = np.einsum("iij->ij", gv_operator).real
deg_sets = degenerate_sets(self._frequencies[irgp])
deg_sets = degenerate_sets(frequencies[irgp])
# group velocities in the degenerate subspace are obtained diagonalizing the
# velocity operator in the subspace of degeneracy.
for id_dir in range(3):
@ -141,14 +184,14 @@ class ConductivityWignerMixIn(HeatCapacityMixIn):
# Sum all vxv at k*
for j, vxv in enumerate(([0, 0], [1, 1], [2, 2], [1, 2], [0, 2], [0, 1])):
# self._gv_sum2[i_data, :, j] = gv_by_gv_tensor[:, vxv[0], vxv[1]]
# self._gv_by_gv[i_data, :, j] = gv_by_gv_tensor[:, vxv[0], vxv[1]]
# here it is storing the 6 independent components of the v^i x v^j tensor
# i_data is the q-point index
# j indexes the 6 independent component of the symmetric tensor v^i x v^j
self._gv_operator_sum2[i_data, :, :, j] = gv_by_gv_operator_tensor[
:, :, vxv[0], vxv[1]
]
# self._gv_sum2[i_data, :, j] = gv_by_gv_tensor[:, vxv[0], vxv[1]]
# self._gv_by_gv[i_data, :, j] = gv_by_gv_tensor[:, vxv[0], vxv[1]]
def _get_gv_by_gv_operator(self, i_irgp, i_data):
if self._is_kappa_star:
@ -211,13 +254,23 @@ class ConductivityWignerMixIn(HeatCapacityMixIn):
return gv_by_gv_operator, order_kstar
def _allocate_values(self):
super()._allocate_values()
def get_conversion_factor_WTE(volume):
"""Return conversion factor of thermal conductivity."""
return (
(get_physical_units().THz * get_physical_units().Angstrom)
** 2 # --> group velocity
* get_physical_units().EV # --> specific heat is in eV/
* get_physical_units().Hbar # --> transform lorentzian_div_hbar from eV^-1 to s
/ (volume * get_physical_units().Angstrom ** 3)
) # --> unit cell volume
num_band0 = len(self._pp.band_indices)
if self._is_reducible_collision_matrix:
num_grid_points = np.prod(self._pp.mesh_numbers)
else:
num_grid_points = len(self._grid_points)
num_band = len(self._pp.primitive) * 3
self._gv_operator = np.zeros(
(num_grid_points, num_band0, num_band, 3),
order="C",
dtype=self._complex_dtype,
)
self._gv_operator_sum2 = np.zeros(
(num_grid_points, num_band0, num_band, 6),
order="C",
dtype=self._complex_dtype,
)

64
phono3py/conductivity/wigner_direct_solution.py
View File

@ -42,13 +42,13 @@ from phono3py.conductivity.direct_solution import (
diagonalize_collision_matrix,
)
from phono3py.conductivity.wigner_base import (
ConductivityWignerMixIn,
ConductivityWignerComponents,
get_conversion_factor_WTE,
)
from phono3py.phonon3.interaction import Interaction
class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
class ConductivityWignerLBTE(ConductivityLBTEBase):
"""Class of Wigner lattice thermal conductivity under direct-solution.
Authors
@ -104,7 +104,6 @@ class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
solve_collective_phonon=solve_collective_phonon,
is_reducible_collision_matrix=is_reducible_collision_matrix,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_full_pp=is_full_pp,
read_pp=read_pp,
pp_filename=pp_filename,
@ -114,10 +113,49 @@ class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
log_level=log_level,
lang=lang,
)
self._conversion_factor_WTE = get_conversion_factor_WTE(
self._pp.primitive.volume
)
self._conductivity_components = ConductivityWignerComponents(
self._pp,
self._grid_points,
self._grid_weights,
self._point_operations,
self._rotations_cartesian,
temperatures=self._temperatures,
is_kappa_star=self._is_kappa_star,
gv_delta_q=gv_delta_q,
is_reducible_collision_matrix=self._is_reducible_collision_matrix,
log_level=self._log_level,
)
@property
def kappa_TOT_RTA(self):
"""Return kappa."""
return self._kappa_TOT_RTA
@property
def kappa_P_RTA(self):
"""Return kappa."""
return self._kappa_P_RTA
@property
def kappa_C(self):
"""Return kappa."""
return self._kappa_C
@property
def mode_kappa_P_RTA(self):
"""Return mode_kappa."""
return self._mode_kappa_P_RTA
@property
def mode_kappa_C(self):
"""Return mode_kappa."""
return self._mode_kappa_C
@property
def kappa_TOT_exact(self):
"""Return kappa."""
@ -133,6 +171,14 @@ class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
"""Return mode_kappa."""
return self._mode_kappa_P_exact
def _set_cv(self, i_gp, i_data):
"""Set cv for conductivity components."""
self._conductivity_components.set_heat_capacities(i_gp, i_data)
def _set_velocities(self, i_gp, i_data):
"""Set velocities for conductivity components."""
self._conductivity_components.set_velocities(i_gp, i_data)
def _allocate_local_values(self, num_grid_points):
"""Allocate grid point local arrays."""
num_band0 = len(self._pp.band_indices)
@ -163,12 +209,6 @@ class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
)
complex_dtype = "c%d" % (np.dtype("double").itemsize * 2)
self._gv_operator = np.zeros(
(num_grid_points, num_band0, nat3, 3), order="C", dtype=complex_dtype
)
self._gv_operator_sum2 = np.zeros(
(num_grid_points, num_band0, nat3, 6), order="C", dtype=complex_dtype
)
# one more index because we have off-diagonal terms (second index not
# parallelized)
self._mode_kappa_C = np.zeros(
@ -279,7 +319,7 @@ class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
def _set_kappa_C_ir_colmat(self, i_sigma, i_temp):
"""Calculate coherence term of the thermal conductivity."""
N = self._num_sampling_grid_points
N = self.number_of_sampling_grid_points
num_band = len(self._pp.primitive) * 3
# num_ir_grid_points = len(self._ir_grid_points)
THzToEv = get_physical_units().THzToEv
@ -287,7 +327,7 @@ class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
# linewidths at qpoint i, sigma i_sigma, and temperature i_temp
g = self._get_main_diagonal(i, i_sigma, i_temp) * 2.0 # linewidth (FWHM)
frequencies = self._frequencies[gp]
cv = self._cv[i_temp, i, :]
cv = self._conductivity_components.mode_heat_capacities[i_temp, i, :]
for s1 in range(num_band):
for s2 in range(num_band):
hbar_omega_eV_s1 = (
@ -317,7 +357,7 @@ class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
* (cv[s1] / hbar_omega_eV_s1 + cv[s2] / hbar_omega_eV_s2)
)
self._mode_kappa_C[i_sigma, i_temp, i, s1, s2] = (
(self._gv_operator_sum2[i, s1, s2])
(self._conductivity_components.gv_by_gv_operator[i, s1, s2])
* prefactor
* lorentzian_divided_by_hbar
* self._conversion_factor_WTE

112
phono3py/conductivity/wigner_rta.py
View File

@ -34,18 +34,20 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import numpy as np
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.rta_base import ConductivityRTABase
from phono3py.conductivity.wigner_base import (
ConductivityWignerMixIn,
ConductivityWignerComponents,
get_conversion_factor_WTE,
)
from phono3py.phonon3.interaction import Interaction
class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
class ConductivityWignerRTA(ConductivityRTABase):
"""Class of Wigner lattice thermal conductivity under RTA.
Authors
@ -57,24 +59,24 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
def __init__(
self,
interaction: Interaction,
grid_points=None,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
is_isotope=False,
mass_variances=None,
boundary_mfp=None, # in micrometer
use_ave_pp=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
read_pp=False,
store_pp=False,
pp_filename=None,
is_N_U=False,
is_gamma_detail=False,
is_frequency_shift_by_bubble=False,
log_level=0,
grid_points: np.ndarray | None = None,
temperatures: list | np.ndarray | None = None,
sigmas: list | np.ndarray | None = None,
sigma_cutoff: float | None = None,
is_isotope: bool = False,
mass_variances: list | np.ndarray | None = None,
boundary_mfp: float | None = None, # in micrometer
use_ave_pp: bool = False,
is_kappa_star: bool = True,
gv_delta_q: float | None = None,
is_full_pp: bool = False,
read_pp: bool = False,
store_pp: bool = False,
pp_filename: float | None = None,
is_N_U: bool = False,
is_gamma_detail: bool = False,
is_frequency_shift_by_bubble: bool = False,
log_level: int = 0,
):
"""Init method."""
self._cv = None
@ -98,7 +100,6 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
boundary_mfp=boundary_mfp,
use_ave_pp=use_ave_pp,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_full_pp=is_full_pp,
read_pp=read_pp,
store_pp=store_pp,
@ -112,6 +113,43 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
self._pp.primitive.volume
)
self._conductivity_components = ConductivityWignerComponents(
self._pp,
self._grid_points,
self._grid_weights,
self._point_operations,
self._rotations_cartesian,
temperatures=self._temperatures,
is_kappa_star=self._is_kappa_star,
gv_delta_q=gv_delta_q,
log_level=self._log_level,
)
@property
def kappa_TOT_RTA(self):
"""Return kappa."""
return self._kappa_TOT_RTA
@property
def kappa_P_RTA(self):
"""Return kappa."""
return self._kappa_P_RTA
@property
def kappa_C(self):
"""Return kappa."""
return self._kappa_C
@property
def mode_kappa_P_RTA(self):
"""Return mode_kappa."""
return self._mode_kappa_P_RTA
@property
def mode_kappa_C(self):
"""Return mode_kappa."""
return self._mode_kappa_C
def set_kappa_at_sigmas(self):
"""Calculate the Wigner thermal conductivity.
@ -120,8 +158,10 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
"""
num_band = len(self._pp.primitive) * 3
THzToEv = get_physical_units().THzToEv
gv_by_gv = self._conductivity_components.gv_by_gv_operator
cv = self._conductivity_components.mode_heat_capacities
for i, _ in enumerate(self._grid_points):
cv = self._cv[:, i, :]
gp = self._grid_points[i]
frequencies = self._frequencies[gp]
# Kappa
@ -156,15 +196,15 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
0.25
* (hbar_omega_eV_s1 + hbar_omega_eV_s2)
* (
cv[k, s1] / hbar_omega_eV_s1
+ cv[k, s2] / hbar_omega_eV_s2
cv[k, i, s1] / hbar_omega_eV_s1
+ cv[k, i, s2] / hbar_omega_eV_s2
)
)
if np.abs(frequencies[s1] - frequencies[s2]) < 1e-4:
# degenerate or diagonal s1=s2 modes contribution
# determine k_P
contribution = (
(self._gv_operator_sum2[i, s1, s2])
(gv_by_gv[i, s1, s2])
* prefactor
* lorentzian_divided_by_hbar
* self._conversion_factor_WTE
@ -182,7 +222,7 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
# conductivity
else:
self._mode_kappa_C[j, k, i, s1, s2] += (
(self._gv_operator_sum2[i, s1, s2])
(gv_by_gv[i, s1, s2])
* prefactor
* lorentzian_divided_by_hbar
* self._conversion_factor_WTE
@ -191,13 +231,21 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
elif s1 == s2:
self._num_ignored_phonon_modes[j, k] += 1
N = self._num_sampling_grid_points
N = self.number_of_sampling_grid_points
self._kappa_P_RTA = self._mode_kappa_P_RTA.sum(axis=2).sum(axis=2) / N
#
self._kappa_C = self._mode_kappa_C.sum(axis=2).sum(axis=2).sum(axis=2) / N
#
self._kappa_TOT_RTA = self._kappa_P_RTA + self._kappa_C
def _set_cv(self, i_gp, i_data):
"""Set cv for conductivity components."""
self._conductivity_components.set_heat_capacities(i_gp, i_data)
def _set_velocities(self, i_gp, i_data):
"""Set velocities for conductivity components."""
self._conductivity_components.set_velocities(i_gp, i_data)
def _allocate_values(self):
super()._allocate_values()
@ -206,16 +254,6 @@ class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
num_temp = len(self._temperatures)
nat3 = len(self._pp.primitive) * 3
self._cv = np.zeros(
(num_temp, num_grid_points, num_band0), order="C", dtype="double"
)
self._gv_operator = np.zeros(
(num_grid_points, num_band0, nat3, 3), order="C", dtype=self._complex_dtype
)
self._gv_operator_sum2 = np.zeros(
(num_grid_points, num_band0, nat3, 6), order="C", dtype=self._complex_dtype
)
# kappa* and mode_kappa* are accessed when all bands exist, i.e.,
# num_band0==num_band.
self._kappa_TOT_RTA = np.zeros(

115
phono3py/cui/create_force_constants.py
View File

@ -51,12 +51,9 @@ from phonopy.harmonic.force_constants import (
symmetrize_force_constants,
)
from phonopy.interface.calculator import get_calculator_physical_units
from phonopy.interface.fc_calculator import fc_calculator_names
from phonopy.interface.pypolymlp import PypolymlpParams, parse_mlp_params
from phonopy.interface.symfc import parse_symfc_options
from phono3py import Phono3py
from phono3py.cui.settings import Phono3pySettings
from phono3py.cui.show_log import show_phono3py_force_constants_settings
from phono3py.file_IO import (
get_length_of_first_line,
@ -67,7 +64,11 @@ from phono3py.file_IO import (
write_fc2_to_hdf5,
write_fc3_to_hdf5,
)
from phono3py.interface.fc_calculator import extract_fc2_fc3_calculators
from phono3py.interface.fc_calculator import (
determine_cutoff_pair_distance,
extract_fc2_fc3_calculators,
get_fc_calculator_params,
)
from phono3py.interface.phono3py_yaml import Phono3pyYaml
from phono3py.phonon3.dataset import forces_in_dataset
from phono3py.phonon3.fc3 import (
@ -77,22 +78,6 @@ from phono3py.phonon3.fc3 import (
)
def get_cutoff_pair_distance(settings: Phono3pySettings) -> Optional[float]:
"""Return cutoff_pair_distance from settings."""
_, fc_calculator_options = get_fc_calculator_params(settings)
if settings.cutoff_pair_distance is None:
cutoff = parse_symfc_options(
extract_fc2_fc3_calculators(fc_calculator_options, 3)
).get("cutoff")
if cutoff is None:
cutoff_pair_distance = None
else:
cutoff_pair_distance = cutoff.get(3)
else:
cutoff_pair_distance = settings.cutoff_pair_distance
return cutoff_pair_distance
def create_phono3py_force_constants(
phono3py: Phono3py,
settings,
@ -131,7 +116,10 @@ def create_phono3py_force_constants(
pass
else:
(fc_calculator, fc_calculator_options) = get_fc_calculator_params(
settings, log_level=(not settings.read_fc3) * 1
settings.fc_calculator,
settings.fc_calculator_options,
settings.cutoff_pair_distance,
log_level=(not settings.read_fc3) * 1,
)
if settings.read_fc3:
_read_phono3py_fc3(phono3py, symmetrize_fc3r, input_filename, log_level)
@ -154,6 +142,7 @@ def create_phono3py_force_constants(
fc_calculator_options, 3
),
)
assert phono3py.fc3 is not None, "fc3 is not set."
cutoff_distance = settings.cutoff_fc3_distance
if cutoff_distance is not None and cutoff_distance > 0:
@ -169,11 +158,13 @@ def create_phono3py_force_constants(
else:
filename = "fc3." + output_filename + ".hdf5"
if log_level:
print('Writing fc3 to "%s".' % filename)
print(f'Writing fc3 to "{filename}".')
write_fc3_to_hdf5(
phono3py.fc3,
fc3_nonzero_indices=phono3py.fc3_nonzero_indices,
filename=filename,
p2s_map=phono3py.primitive.p2s_map,
fc3_cutoff=phono3py.fc3_cutoff,
compression=settings.hdf5_compression,
)
@ -307,64 +298,6 @@ def parse_forces(
return dataset
def get_fc_calculator_params(
settings: Phono3pySettings, log_level: int = 0
) -> tuple[Optional[str], Optional[str]]:
"""Return fc_calculator and fc_calculator_params from settings."""
fc_calculator = None
fc_calculator_list = []
if settings.fc_calculator is not None:
for fc_calculatr_str in settings.fc_calculator.split("|"):
if fc_calculatr_str == "": # No external calculator
fc_calculator_list.append(fc_calculatr_str.lower())
elif fc_calculatr_str.lower() in fc_calculator_names:
fc_calculator_list.append(fc_calculatr_str.lower())
if fc_calculator_list:
fc_calculator = "|".join(fc_calculator_list)
fc_calculator_options = settings.fc_calculator_options
if settings.cutoff_pair_distance:
if fc_calculator_list and fc_calculator_list[-1] in ("alm", "symfc"):
if fc_calculator_list[-1] == "alm":
cutoff_str = f"-1 {settings.cutoff_pair_distance}"
if fc_calculator_list[-1] == "symfc":
cutoff_str = f"{settings.cutoff_pair_distance}"
fc_calculator_options = _set_cutoff_in_fc_calculator_options(
fc_calculator_options,
cutoff_str,
log_level,
)
return fc_calculator, fc_calculator_options
def _set_cutoff_in_fc_calculator_options(
fc_calculator_options: Optional[str],
cutoff_str: str,
log_level: int,
):
str_appended = f"cutoff={cutoff_str}"
calc_opts = fc_calculator_options
if calc_opts is None:
calc_opts = "|"
if "|" in calc_opts:
calc_opts_fc2, calc_opts_fc3 = [v.strip() for v in calc_opts.split("|")][:2]
else:
calc_opts_fc2 = calc_opts
calc_opts_fc3 = calc_opts
if calc_opts_fc3 == "":
calc_opts_fc3 += f"{str_appended}"
if log_level:
print(f'Set "{str_appended}" to fc_calculator_options for fc3.')
elif "cutoff" not in calc_opts_fc3:
calc_opts_fc3 += f", {str_appended}"
if log_level:
print(f'Appended "{str_appended}" to fc_calculator_options for fc3.')
return f"{calc_opts_fc2}|{calc_opts_fc3}"
def _read_phono3py_fc3(phono3py: Phono3py, symmetrize_fc3r, input_filename, log_level):
if input_filename is None:
filename = "fc3.hdf5"
@ -524,13 +457,19 @@ def run_pypolymlp_to_compute_forces(
number_of_snapshots: Optional[int] = None,
random_seed: Optional[int] = None,
prepare_dataset: bool = False,
fc_calculator: Optional[str] = None,
fc_calculator_options: Optional[str] = None,
cutoff_pair_distance: Optional[float] = None,
symfc_memory_size: Optional[float] = None,
mlp_filename: Optional[str] = None,
log_level: int = 0,
):
"""Run pypolymlp to compute forces."""
if log_level:
import pypolymlp
print("-" * 29 + " pypolymlp start " + "-" * 30)
print("Pypolymlp version", pypolymlp.__version__)
print("Pypolymlp is a generator of polynomial machine learning potentials.")
print("Please cite the paper: A. Seko, J. Appl. Phys. 133, 011101 (2023).")
print("Pypolymlp is developed at https://github.com/sekocha/pypolymlp.")
@ -540,7 +479,6 @@ def run_pypolymlp_to_compute_forces(
_mlp_filename_list = list(pathlib.Path().glob(f"{mlp_filename}*"))
if _mlp_filename_list:
_mlp_filename = _mlp_filename_list[0]
print(_mlp_filename)
if _mlp_filename.suffix not in [
".yaml",
".pmlp",
@ -601,6 +539,18 @@ def run_pypolymlp_to_compute_forces(
"0"
).rstrip(".")
)
cutoff_pair_distance = determine_cutoff_pair_distance(
fc_calculator=fc_calculator,
fc_calculator_options=fc_calculator_options,
cutoff_pair_distance=cutoff_pair_distance,
symfc_memory_size=symfc_memory_size,
random_displacements=number_of_snapshots,
supercell=ph3py.supercell,
primitive=ph3py.primitive,
symmetry=ph3py.symmetry,
log_level=log_level,
)
ph3py.generate_displacements(
distance=_displacement_distance,
cutoff_pair_distance=cutoff_pair_distance,
@ -633,7 +583,10 @@ def run_pypolymlp_to_compute_phonon_forces(
"""Run pypolymlp to compute phonon forces."""
if ph3py.phonon_mlp_dataset is not None:
if log_level:
import pypolymlp
print("-" * 29 + " pypolymlp start " + "-" * 30)
print(f"Pypolymlp version {pypolymlp.__version__}")
print("Pypolymlp is a generator of polynomial machine learning potentials.")
print("Please cite the paper: A. Seko, J. Appl. Phys. 133, 011101 (2023).")
print("Pypolymlp is developed at https://github.com/sekocha/pypolymlp.")

87
phono3py/cui/create_supercells.py
View File

@ -45,6 +45,7 @@ from phono3py.interface.calculator import (
get_additional_info_to_write_supercells,
get_default_displacement_distance,
)
from phono3py.interface.fc_calculator import determine_cutoff_pair_distance
def create_phono3py_supercells(
@ -57,7 +58,7 @@ def create_phono3py_supercells(
"""Create displacements and supercells.
Distance unit used is that for the calculator interface.
The default unit is Angstron.
The default unit is Angstrom.
"""
optional_structure_info = cell_info["optional_structure_info"]
@ -66,7 +67,7 @@ def create_phono3py_supercells(
distance = get_default_displacement_distance(interface_mode)
else:
distance = settings.displacement_distance
phono3py = Phono3py(
ph3 = Phono3py(
cell_info["unitcell"],
cell_info["supercell_matrix"],
primitive_matrix=cell_info["primitive_matrix"],
@ -75,9 +76,34 @@ def create_phono3py_supercells(
symprec=symprec,
calculator=interface_mode,
)
phono3py.generate_displacements(
distance=distance,
if log_level:
print("")
print('Unit cell was read from "%s".' % optional_structure_info[0])
print("-" * 32 + " unit cell " + "-" * 33) # 32 + 11 + 33 = 76
print_cell(ph3.unitcell)
print("-" * 76)
print_supercell_matrix(ph3.supercell_matrix, ph3.phonon_supercell_matrix)
if ph3.primitive_matrix is not None:
print("Primitive matrix:")
for v in ph3.primitive_matrix:
print(" %s" % v)
print("Displacement distance: %s" % distance)
cutoff_pair_distance = determine_cutoff_pair_distance(
fc_calculator=settings.fc_calculator,
fc_calculator_options=settings.fc_calculator_options,
cutoff_pair_distance=settings.cutoff_pair_distance,
symfc_memory_size=settings.symfc_memory_size,
random_displacements=settings.random_displacements,
supercell=ph3.supercell,
primitive=ph3.primitive,
symmetry=ph3.symmetry,
log_level=log_level,
)
ph3.generate_displacements(
distance=distance,
cutoff_pair_distance=cutoff_pair_distance,
is_plusminus=settings.is_plusminus_displacement,
is_diagonal=settings.is_diagonal_displacement,
number_of_snapshots=settings.random_displacements,
@ -88,44 +114,29 @@ def create_phono3py_supercells(
settings.random_displacements_fc2
or settings.phonon_supercell_matrix is not None
):
phono3py.generate_fc2_displacements(
ph3.generate_fc2_displacements(
distance=distance,
is_plusminus=settings.is_plusminus_displacement_fc2,
number_of_snapshots=settings.random_displacements_fc2,
random_seed=settings.random_seed,
)
if log_level:
print("")
print('Unit cell was read from "%s".' % optional_structure_info[0])
print("-" * 32 + " unit cell " + "-" * 33) # 32 + 11 + 33 = 76
print_cell(phono3py.unitcell)
print("-" * 76)
print_supercell_matrix(
phono3py.supercell_matrix, phono3py.phonon_supercell_matrix
)
if phono3py.primitive_matrix is not None:
print("Primitive matrix:")
for v in phono3py.primitive_matrix:
print(" %s" % v)
print("Displacement distance: %s" % distance)
ids = []
disp_cells = []
for i, cell in enumerate(phono3py.supercells_with_displacements):
for i, cell in enumerate(ph3.supercells_with_displacements):
if cell is not None:
ids.append(i + 1)
disp_cells.append(cell)
additional_info = get_additional_info_to_write_supercells(
interface_mode, phono3py.supercell_matrix
interface_mode, ph3.supercell_matrix
)
additional_info["supercell_matrix"] = phono3py.supercell_matrix
additional_info["supercell_matrix"] = ph3.supercell_matrix
write_supercells_with_displacements(
interface_mode,
phono3py.supercell,
ph3.supercell,
disp_cells,
optional_structure_info,
displacement_ids=ids,
@ -134,24 +145,22 @@ def create_phono3py_supercells(
)
if log_level:
num_disps = len(phono3py.supercells_with_displacements)
num_disps = len(ph3.supercells_with_displacements)
num_disp_files = len(disp_cells)
print("Number of displacements: %d" % num_disps)
if settings.cutoff_pair_distance is not None:
print(
"Cutoff distance for displacements: %s" % settings.cutoff_pair_distance
)
print("Number of displacement supercell files created: %d" % num_disp_files)
print(f"Number of displacements: {num_disps}")
if cutoff_pair_distance is not None:
print(f"Cutoff distance for displacements: {cutoff_pair_distance}")
print(f"Number of displacement supercell files created: {num_disp_files}")
if phono3py.phonon_supercell_matrix is not None:
num_disps = len(phono3py.phonon_supercells_with_displacements)
if ph3.phonon_supercell_matrix is not None:
num_disps = len(ph3.phonon_supercells_with_displacements)
additional_info = get_additional_info_to_write_fc2_supercells(
interface_mode, phono3py.phonon_supercell_matrix
interface_mode, ph3.phonon_supercell_matrix
)
write_supercells_with_displacements(
interface_mode,
phono3py.phonon_supercell,
phono3py.phonon_supercells_with_displacements,
ph3.phonon_supercell,
ph3.phonon_supercells_with_displacements,
optional_structure_info,
zfill_width=5,
additional_info=additional_info,
@ -165,11 +174,11 @@ def create_phono3py_supercells(
n_pure_trans = sum(
[
(r == identity).all()
for r in phono3py.symmetry.symmetry_operations["rotations"]
for r in ph3.symmetry.symmetry_operations["rotations"]
]
)
if len(phono3py.supercell) // len(phono3py.primitive) != n_pure_trans:
if len(ph3.supercell) // len(ph3.primitive) != n_pure_trans:
print("*" * 72)
print(
"Note: "
@ -177,4 +186,4 @@ def create_phono3py_supercells(
)
print("*" * 72)
return phono3py
return ph3

53
phono3py/cui/load.py
View File

@ -42,6 +42,7 @@ from typing import Optional, Union
import numpy as np
import phonopy.cui.load_helper as load_helper
from numpy.typing import NDArray
from phonopy.harmonic.force_constants import show_drift_force_constants
from phonopy.interface.calculator import get_calculator_physical_units
from phonopy.physical_units import get_physical_units
@ -392,12 +393,10 @@ def load_fc2_and_fc3(
):
"""Set force constants."""
if fc3_filename is not None or pathlib.Path("fc3.hdf5").exists():
fc3 = _load_fc3(ph3py, fc3_filename=fc3_filename, log_level=log_level)
ph3py.fc3 = fc3
_load_fc3(ph3py, fc3_filename=fc3_filename, log_level=log_level)
if fc2_filename is not None or pathlib.Path("fc2.hdf5").exists():
fc2 = _load_fc2(ph3py, fc2_filename=fc2_filename, log_level=log_level)
ph3py.fc2 = fc2
_load_fc2(ph3py, fc2_filename=fc2_filename, log_level=log_level)
def load_dataset_and_phonon_dataset(
@ -493,30 +492,40 @@ def compute_force_constants_from_datasets(
def _load_fc3(
ph3py: Phono3py,
fc3_filename: Optional[os.PathLike] = None,
fc3_filename: str | os.PathLike | None = None,
log_level: int = 0,
) -> np.ndarray:
):
p2s_map = ph3py.primitive.p2s_map
if fc3_filename is None:
_fc3_filename = "fc3.hdf5"
else:
_fc3_filename = fc3_filename
fc3 = read_fc3_from_hdf5(filename=_fc3_filename, p2s_map=p2s_map)
_check_fc3_shape(ph3py, fc3, filename=_fc3_filename)
if log_level:
print(f'fc3 was read from "{_fc3_filename}".')
return fc3
if isinstance(fc3, dict):
# fc3 is read from a file with type-1 format.
assert "fc3" in fc3
_check_fc3_shape(ph3py, fc3["fc3"], filename=_fc3_filename)
ph3py.fc3 = fc3["fc3"]
assert "fc3_nonzero_indices" in fc3
ph3py.fc3_nonzero_indices = fc3["fc3_nonzero_indices"]
if log_level:
print(f'fc3 and fc3 nonzero indices were read from "{_fc3_filename}".')
else:
_check_fc3_shape(ph3py, fc3, filename=_fc3_filename)
ph3py.fc3 = fc3
if log_level:
print(f'fc3 was read from "{_fc3_filename}".')
def _select_and_load_dataset(
ph3py: Phono3py,
ph3py_yaml: Optional[Phono3pyYaml] = None,
forces_fc3_filename: Optional[Union[os.PathLike, Sequence]] = None,
phono3py_yaml_filename: Optional[os.PathLike] = None,
cutoff_pair_distance: Optional[float] = None,
calculator: Optional[str] = None,
ph3py_yaml: Phono3pyYaml | None = None,
forces_fc3_filename: os.PathLike | Sequence | None = None,
phono3py_yaml_filename: os.PathLike | None = None,
cutoff_pair_distance: float | None = None,
calculator: str | None = None,
log_level: int = 0,
) -> Optional[dict]:
) -> dict | None:
dataset = None
if (
ph3py_yaml is not None
@ -563,8 +572,8 @@ def _select_and_load_dataset(
def _load_fc2(
ph3py: Phono3py, fc2_filename: Optional[os.PathLike] = None, log_level: int = 0
) -> np.ndarray:
ph3py: Phono3py, fc2_filename: os.PathLike | None = None, log_level: int = 0
):
phonon_p2s_map = ph3py.phonon_primitive.p2s_map
if fc2_filename is None:
_fc2_filename = "fc2.hdf5"
@ -574,7 +583,7 @@ def _load_fc2(
_check_fc2_shape(ph3py, fc2, filename=_fc2_filename)
if log_level:
print(f'fc2 was read from "{_fc2_filename}".')
return fc2
ph3py.fc2 = fc2
def _select_and_load_phonon_dataset(
@ -669,7 +678,7 @@ def _get_dataset_for_fc2(
return dataset
def _check_fc2_shape(ph3py: Phono3py, fc2, filename="fc2.hdf5"):
def _check_fc2_shape(ph3py: Phono3py, fc2, filename: str | os.PathLike = "fc2.hdf5"):
if ph3py.phonon_supercell_matrix is None:
smat = ph3py.supercell_matrix
else:
@ -677,7 +686,9 @@ def _check_fc2_shape(ph3py: Phono3py, fc2, filename="fc2.hdf5"):
_check_fc_shape(ph3py, fc2, smat, filename)
def _check_fc3_shape(ph3py: Phono3py, fc3, filename="fc3.hdf5"):
def _check_fc3_shape(
ph3py: Phono3py, fc3: NDArray, filename: str | os.PathLike = "fc3.hdf5"
):
smat = ph3py.supercell_matrix
_check_fc_shape(ph3py, fc3, smat, filename)

21
phono3py/cui/phono3py_argparse.py
View File

@ -288,13 +288,6 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
default=False,
help="Read third order force constants",
)
parser.add_argument(
"--v2",
dest="is_fc3_r0_average",
action="store_false",
default=True,
help="Take average in fc3-r2q transformation around three atoms",
)
parser.add_argument(
"--fc-calc",
"--fc-calculator",
@ -778,6 +771,13 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
default=None,
help="Show symfc memory usage with respect to cutoff distance",
)
parser.add_argument(
"--symfc-memsize",
dest="symfc_memory_size",
type=float,
default=None,
help="Memory size to estimate cutoff pair distance (GB)",
)
parser.add_argument(
"--spf",
dest="is_spectral_function",
@ -862,6 +862,13 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
default=False,
help="Detailed run-time information is displayed",
)
parser.add_argument(
"--v2",
dest="is_fc3_r0_average",
action="store_false",
default=True,
help="Take average in fc3-r2q transformation around three atoms",
)
parser.add_argument(
"--wgp",
"--write-grid-points",

54
phono3py/cui/phono3py_script.py
View File

@ -59,6 +59,7 @@ from phonopy.exception import ForceCalculatorRequiredError
from phonopy.file_IO import is_file_phonopy_yaml
from phonopy.harmonic.force_constants import show_drift_force_constants
from phonopy.interface.calculator import get_calculator_physical_units
from phonopy.interface.symfc import estimate_symfc_cutoff_from_memsize
from phonopy.phonon.band_structure import get_band_qpoints
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import isclose as cells_isclose
@ -66,7 +67,6 @@ from phonopy.structure.cells import isclose as cells_isclose
from phono3py import Phono3py, Phono3pyIsotope, Phono3pyJointDos
from phono3py.cui.create_force_constants import (
create_phono3py_force_constants,
get_cutoff_pair_distance,
get_fc_calculator_params,
run_pypolymlp_to_compute_forces,
)
@ -95,7 +95,7 @@ from phono3py.file_IO import (
write_fc3_to_hdf5,
write_phonon_to_hdf5,
)
from phono3py.interface.fc_calculator import estimate_symfc_memory_usage
from phono3py.interface.fc_calculator import determine_cutoff_pair_distance
from phono3py.interface.phono3py_yaml import Phono3pyYaml
from phono3py.phonon.grid import get_grid_point_from_address, get_ir_grid_points
from phono3py.phonon3.dataset import forces_in_dataset
@ -593,14 +593,33 @@ def _store_force_constants(ph3py: Phono3py, settings: Phono3pySettings, log_leve
if log_level:
print("-" * 29 + " Force constants " + "-" * 30)
load_fc2_and_fc3(ph3py, log_level=log_level)
read_fc3 = ph3py.fc3 is not None
read_fc2 = ph3py.fc2 is not None
load_fc2_and_fc3(ph3py, log_level=log_level)
cutoff_pair_distance = None
if settings.use_pypolymlp:
cutoff_pair_distance = ph3py.dataset.get("cutoff_distance")
if cutoff_pair_distance is None:
cutoff_pair_distance = determine_cutoff_pair_distance(
fc_calculator=settings.fc_calculator,
fc_calculator_options=settings.fc_calculator_options,
cutoff_pair_distance=settings.cutoff_pair_distance,
symfc_memory_size=settings.symfc_memory_size,
random_displacements=settings.random_displacements,
supercell=ph3py.supercell,
primitive=ph3py.primitive,
log_level=log_level,
)
if cutoff_pair_distance is None and ph3py.dataset is not None:
cutoff_pair_distance = ph3py.dataset.get("cutoff_distance")
cutoff_pair_distance = get_cutoff_pair_distance(settings)
(fc_calculator, fc_calculator_options) = get_fc_calculator_params(
settings, log_level=(not read_fc3) * 1
settings.fc_calculator,
settings.fc_calculator_options,
cutoff_pair_distance,
log_level=(not read_fc3) * 1,
)
try:
compute_force_constants_from_datasets(
@ -660,7 +679,9 @@ def _store_force_constants(ph3py: Phono3py, settings: Phono3pySettings, log_leve
if not read_fc3:
write_fc3_to_hdf5(
ph3py.fc3,
fc3_nonzero_indices=ph3py.fc3_nonzero_indices,
p2s_map=ph3py.primitive.p2s_map,
fc3_cutoff=ph3py.fc3_cutoff,
compression=settings.hdf5_compression,
)
if log_level:
@ -668,7 +689,7 @@ def _store_force_constants(ph3py: Phono3py, settings: Phono3pySettings, log_leve
if not read_fc2:
write_fc2_to_hdf5(
ph3py.fc2,
p2s_map=ph3py.primitive.p2s_map,
p2s_map=ph3py.phonon_primitive.p2s_map,
physical_unit="eV/angstrom^2",
compression=settings.hdf5_compression,
)
@ -1048,20 +1069,11 @@ def main(**argparse_control):
###############################
if settings.show_symfc_memory_usage and load_phono3py_yaml:
print("Quick estimation of memory size required for solving fc3 by symfc")
vecs, _ = ph3py.primitive.get_smallest_vectors()
dists = np.unique(
np.round(np.linalg.norm(vecs @ ph3py.primitive.cell, axis=-1), decimals=1)
)
print("cutoff memsize")
print("------ -------")
for cutoff in dists[1:] + 0.1:
memsize, memsize2 = estimate_symfc_memory_usage(
ph3py.supercell, ph3py.symmetry, cutoff
)
print(
f"{cutoff:5.1f} {memsize + memsize2:6.2f} GB "
f"({memsize:.2f}+{memsize2:.2f})"
)
estimate_symfc_cutoff_from_memsize(
ph3py.supercell, ph3py.primitive, ph3py.symmetry, 3, verbose=True
)
if log_level:
print_end_phono3py()
@ -1179,14 +1191,16 @@ def main(**argparse_control):
prepare_dataset = (
settings.create_displacements or settings.random_displacements is not None
)
cutoff_pair_distance = get_cutoff_pair_distance(settings)
run_pypolymlp_to_compute_forces(
ph3py,
mlp_params=settings.mlp_params,
displacement_distance=settings.displacement_distance,
number_of_snapshots=settings.random_displacements,
random_seed=settings.random_seed,
cutoff_pair_distance=cutoff_pair_distance,
fc_calculator=settings.fc_calculator,
fc_calculator_options=settings.fc_calculator_options,
cutoff_pair_distance=settings.cutoff_pair_distance,
symfc_memory_size=settings.symfc_memory_size,
prepare_dataset=prepare_dataset,
log_level=log_level,
)

20
phono3py/cui/settings.py
View File

@ -51,6 +51,7 @@ class Phono3pySettings(Settings):
"create_forces_fc3_file": None,
"cutoff_fc3_distance": None,
"cutoff_pair_distance": None,
"emulate_v2": False,
"grid_addresses": None,
"grid_points": None,
"grid_matrix": None,
@ -94,8 +95,8 @@ class Phono3pySettings(Settings):
"scattering_event_class": None, # scattering event class 1 or 2
"sigma_cutoff_width": None,
"solve_collective_phonon": False,
"emulate_v2": False,
"show_symfc_memory_usage": False,
"symfc_memory_size": None,
"subtract_forces": None,
"subtract_forces_fc2": None,
"temperatures": None,
@ -336,6 +337,10 @@ class Phono3pySettings(Settings):
"""Set subtract_forces_fc2."""
self._v["subtract_forces_fc2"] = val
def set_symfc_memory_size(self, val):
"""Set symfc_memory_size."""
self._v["symfc_memory_size"] = val
def set_temperatures(self, val):
"""Set temperatures."""
self._v["temperatures"] = val
@ -612,9 +617,8 @@ class Phono3pyConfParser(ConfParser):
self._confs["scattering_event_class"] = scatt_class
if "sigma_cutoff_width" in self._args:
sigma_cutoff = self._args.sigma_cutoff_width
if sigma_cutoff is not None:
self._confs["sigma_cutoff_width"] = sigma_cutoff
if self._args.sigma_cutoff_width is not None:
self._confs["sigma_cutoff_width"] = self._args.sigma_cutoff_width
if "solve_collective_phonon" in self._args:
if self._args.solve_collective_phonon:
@ -632,6 +636,10 @@ class Phono3pyConfParser(ConfParser):
if self._args.subtract_forces_fc2:
self._confs["subtract_forces_fc2"] = self._args.subtract_forces_fc2
if "symfc_memory_size" in self._args:
if self._args.symfc_memory_size is not None:
self._confs["symfc_memory_size"] = self._args.symfc_memory_size
if "temperatures" in self._args:
if self._args.temperatures is not None:
self._confs["temperatures"] = " ".join(self._args.temperatures)
@ -725,6 +733,7 @@ class Phono3pyConfParser(ConfParser):
"pinv_cutoff",
"pp_conversion_factor",
"sigma_cutoff_width",
"symfc_memory_size",
):
self.set_parameter(conf_key, float(confs[conf_key]))
@ -1067,6 +1076,9 @@ class Phono3pyConfParser(ConfParser):
if "subtract_forces_fc2" in params:
self._settings.set_subtract_forces_fc2(params["subtract_forces_fc2"])
if "symfc_memory_size" in params:
self._settings.set_symfc_memory_size(params["symfc_memory_size"])
# Temperatures for scatterings
if "temperatures" in params:
self._settings.set_temperatures(params["temperatures"])

116
phono3py/file_IO.py
View File

@ -36,12 +36,14 @@
from __future__ import annotations
import os
import pathlib
import warnings
from collections.abc import Sequence
from typing import Optional, TextIO, Union
import h5py
import numpy as np
from numpy.typing import NDArray
from phonopy.cui.load_helper import read_force_constants_from_hdf5
# This import is deactivated for a while.
@ -300,34 +302,54 @@ def _write_FORCES_FC3_typeI(
count += 1
def write_fc3_to_hdf5(fc3, filename="fc3.hdf5", p2s_map=None, compression="gzip"):
def write_fc3_to_hdf5(
fc3: NDArray,
fc3_nonzero_indices: NDArray | None = None,
filename: str = "fc3.hdf5",
p2s_map: NDArray | None = None,
fc3_cutoff: float | None = None,
compression: str = "gzip",
):
"""Write fc3 in fc3.hdf5.
Parameters
----------
force_constants : ndarray
Force constants
shape=(n_satom, n_satom, n_satom, 3, 3, 3) or
(n_patom, n_satom, n_satom,3,3,3), dtype=double
fc3 : ndarray
Force constants shape=(n_satom, n_satom, n_satom, 3, 3, 3) or (n_patom,
n_satom, n_satom,3,3,3), dtype=double
fc3_nonzero_indices : ndarray, optional
Non-zero indices of fc3. shape=(n_satom, n_satom, n_satom) or (n_patom,
n_satom, n_satom), dtype="byte". If this is given, fc3 is in compact
format. Otherwise, it is in full format.
filename : str
Filename to be used.
p2s_map : ndarray, optional
Primitive atom indices in supercell index system
shape=(n_patom,), dtype=intc
Primitive atom indices in supercell index system shape=(n_patom,),
dtype=intc
fc3_cutoff : float, optional
Cutoff distance for fc3.
compression : str or int, optional
h5py's lossless compression filters (e.g., "gzip", "lzf"). None gives
no compression. See the detail at docstring of
h5py.Group.create_dataset. Default is "gzip".
h5py's lossless compression filters (e.g., "gzip", "lzf"). None gives no
compression. See the detail at docstring of h5py.Group.create_dataset.
Default is "gzip".
"""
with h5py.File(filename, "w") as w:
w.create_dataset("version", data=np.bytes_(__version__))
w.create_dataset("fc3", data=fc3, compression=compression)
if fc3_nonzero_indices is not None:
w.create_dataset(
"fc3_nonzero_indices", data=fc3_nonzero_indices, compression=compression
)
if fc3_cutoff is not None:
w.create_dataset("fc3_cutoff", data=fc3_cutoff)
if p2s_map is not None:
w.create_dataset("p2s_map", data=p2s_map)
def read_fc3_from_hdf5(filename="fc3.hdf5", p2s_map=None):
def read_fc3_from_hdf5(
filename: str | os.PathLike = "fc3.hdf5", p2s_map: NDArray | None = None
) -> NDArray | dict:
"""Read fc3 from fc3.hdf5.
fc3 can be in full or compact format. They are distinguished by
@ -339,20 +361,44 @@ def read_fc3_from_hdf5(filename="fc3.hdf5", p2s_map=None):
"""
with h5py.File(filename, "r") as f:
fc3 = f["fc3"][:]
if "fc3" not in f:
raise KeyError(
f"{filename} does not have 'fc3' dataset. "
"This file is not a valid fc3.hdf5."
)
fc3: NDArray = f["fc3"][:] # type: ignore
if fc3.dtype == np.dtype("double") and fc3.flags.c_contiguous:
pass
else:
raise TypeError(
f"{filename} has to be read by h5py as numpy ndarray of "
"dtype='double' and c_contiguous."
)
if "p2s_map" in f:
p2s_map_in_file = f["p2s_map"][:]
check_force_constants_indices(
fc3.shape[:2], p2s_map_in_file, p2s_map, filename
)
if fc3.dtype == np.dtype("double") and fc3.flags.c_contiguous:
fc3_nonzero_indices = None # type: ignore
if "fc3_nonzero_indices" in f:
fc3_nonzero_indices: NDArray = f["fc3_nonzero_indices"][:] # type: ignore
if (
fc3_nonzero_indices.dtype == np.dtype("byte")
and fc3_nonzero_indices.flags.c_contiguous
):
pass
else:
raise TypeError(
f"{filename} has to be read by h5py as numpy ndarray of "
"dtype='byte' and c_contiguous."
)
if fc3_nonzero_indices is None:
return fc3
else:
msg = (
"%s has to be read by h5py as numpy ndarray of "
"dtype='double' and c_contiguous." % filename
)
raise TypeError(msg)
return {"fc3": fc3, "fc3_nonzero_indices": fc3_nonzero_indices}
def write_fc2_to_hdf5(
@ -415,7 +461,7 @@ def read_fc2_from_hdf5(filename="fc2.hdf5", p2s_map=None):
def write_datasets_to_hdf5(
dataset: dict,
phonon_dataset: dict = None,
phonon_dataset: dict | None = None,
filename: str = "datasets.hdf5",
compression: str = "gzip",
):
@ -1090,7 +1136,7 @@ def read_gamma_from_hdf5(
filename=filename,
)
full_filename = "kappa" + suffix + ".hdf5"
if not os.path.exists(full_filename):
if not pathlib.Path(full_filename).exists():
return None, full_filename
read_data = {}
@ -1117,8 +1163,8 @@ def read_collision_from_hdf5(
filename=None,
only_temperatures=False,
verbose=True,
):
"""Read colliction matrix.
) -> tuple:
"""Read collision matrix.
indices : array_like of int
Indices of temperatures.
@ -1137,10 +1183,8 @@ def read_collision_from_hdf5(
filename=filename,
)
full_filename = "collision" + suffix + ".hdf5"
if not os.path.exists(full_filename):
if verbose:
print("%s not found." % full_filename)
return None
if not pathlib.Path(full_filename).exists():
raise FileNotFoundError(f"{full_filename} not found.")
if only_temperatures:
with h5py.File(full_filename, "r") as f:
@ -1290,7 +1334,7 @@ def read_pp_from_hdf5(
filename=None,
verbose=True,
check_consistency=False,
):
) -> tuple:
"""Read ph-ph interaction strength from its hdf5 file."""
suffix = _get_filename_suffix(
mesh,
@ -1300,10 +1344,8 @@ def read_pp_from_hdf5(
filename=filename,
)
full_filename = "pp" + suffix + ".hdf5"
if not os.path.exists(full_filename):
if verbose:
print("%s not found." % full_filename)
return None
if not pathlib.Path(full_filename).exists():
raise FileNotFoundError(f"{full_filename} not found.")
with h5py.File(full_filename, "r") as f:
if "nonzero_pp" in f:
@ -1346,7 +1388,7 @@ def read_pp_from_hdf5(
return pp, g_zero
return None
raise RuntimeError(f"pp data not found in {full_filename}.")
def write_gamma_detail_to_hdf5(
@ -1479,14 +1521,12 @@ def write_phonon_to_hdf5(
return None
def read_phonon_from_hdf5(mesh, filename=None, verbose=True):
def read_phonon_from_hdf5(mesh, filename=None, verbose=True) -> tuple:
"""Read phonon from its hdf5 file."""
suffix = _get_filename_suffix(mesh, filename=filename)
full_filename = "phonon" + suffix + ".hdf5"
if not os.path.exists(full_filename):
if verbose:
print("%s not found." % full_filename)
return None
if not pathlib.Path(full_filename).exists():
raise FileNotFoundError(f"{full_filename} not found.")
with h5py.File(full_filename, "r") as f:
frequencies = np.array(f["frequency"][:], dtype="double", order="C")
@ -1504,8 +1544,6 @@ def read_phonon_from_hdf5(mesh, filename=None, verbose=True):
return frequencies, eigenvectors, grid_address
return None
def write_ir_grid_points(bz_grid, grid_points, grid_weights, primitive_lattice):
"""Write ir-grid-points in yaml."""

322
phono3py/interface/fc_calculator.py
View File

@ -39,8 +39,8 @@ from __future__ import annotations
from typing import Optional, Union
import numpy as np
from phonopy.interface.fc_calculator import FCSolver
from phonopy.interface.symfc import SymfcFCSolver
from phonopy.interface.fc_calculator import FCSolver, fc_calculator_names
from phonopy.interface.symfc import parse_symfc_options, update_symfc_cutoff_by_memsize
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import Primitive
from phonopy.structure.symmetry import Symmetry
@ -49,137 +49,6 @@ from phono3py.phonon3.dataset import get_displacements_and_forces_fc3
from phono3py.phonon3.fc3 import get_fc3
def get_fc3_solver(
supercell: PhonopyAtoms,
primitive: Primitive,
dataset: dict,
fc_calculator: Optional[str] = None,
fc_calculator_options: Optional[str] = None,
is_compact_fc: bool = False,
symmetry: Optional[Symmetry] = None,
log_level: int = 0,
) -> FC3Solver:
"""Return force constants solver for fc3.
Parameters
----------
supercell : PhonopyAtoms
Supercell
primitive : Primitive
Primitive cell
dataset : dict, optional
Dataset that contains displacements, forces, and optionally
energies. Default is None.
fc_calculator : str, optional
Currently only 'alm' is supported. Default is None, meaning invoking
'alm'.
fc_calculator_options : str, optional
This is arbitrary string.
is_compact_fc : bool, optional
If True, force constants are returned in the compact form.
symmetry : Symmetry, optional
Symmetry of supercell. This is used for the traditional and symfc FC
solver. Default is None.
log_level : integer or bool, optional
Verbosity level. False or 0 means quiet. True or 1 means normal level
of log to stdout. 2 gives verbose mode.
Returns
-------
FC3Solver
Force constants solver for fc3.
"""
fc_solver_name = fc_calculator if fc_calculator is not None else "traditional"
fc_solver = FC3Solver(
fc_solver_name,
supercell,
symmetry=symmetry,
dataset=dataset,
is_compact_fc=is_compact_fc,
primitive=primitive,
orders=[2, 3],
options=fc_calculator_options,
log_level=log_level,
)
return fc_solver
def extract_fc2_fc3_calculators(
fc_calculator: Optional[Union[str, dict]],
order: int,
) -> Optional[Union[str, dict]]:
"""Extract fc_calculator and fc_calculator_options for fc2 and fc3.
fc_calculator : str
FC calculator. "|" separates fc2 and fc3. First and last
parts separated correspond to fc2 and fc3 calculators, respectively.
order : int = 2 or 3
2 and 3 indicate fc2 and fc3, respectively.
"""
if isinstance(fc_calculator, dict) or fc_calculator is None:
return fc_calculator
elif isinstance(fc_calculator, str):
if "|" in fc_calculator:
_fc_calculator = fc_calculator.split("|")[order - 2]
if _fc_calculator == "":
return None
return _fc_calculator
else:
if fc_calculator.strip() == "":
return None
return fc_calculator
else:
raise RuntimeError("fc_calculator should be str, dict, or None.")
def update_cutoff_fc_calculator_options(
fc_calc_opts: Optional[Union[str, dict]],
cutoff_pair_distance: Optional[float],
) -> Optional[Union[str, dict]]:
"""Update fc_calculator_options with cutoff distances.
Parameters
----------
fc_calc_opts : str or dict
FC calculator options.
cutoff_pair_distance : float, optional
Cutoff distance for pair interaction.
"""
if cutoff_pair_distance is not None:
if not isinstance(fc_calc_opts, (str, dict)) and fc_calc_opts is not None:
raise RuntimeError("fc_calculator_options should be str, dict, or None.")
if isinstance(fc_calc_opts, dict) and "cutoff" not in fc_calc_opts:
fc_calc_opts["cutoff"] = float(cutoff_pair_distance)
elif isinstance(fc_calc_opts, str) and "cutoff" not in fc_calc_opts:
fc_calc_opts = f"{fc_calc_opts}, cutoff = {cutoff_pair_distance}"
elif fc_calc_opts is None:
fc_calc_opts = f"cutoff = {cutoff_pair_distance}"
return fc_calc_opts
def estimate_symfc_memory_usage(
supercell: PhonopyAtoms, symmetry: Symmetry, cutoff: float, batch_size: int = 100
):
"""Estimate memory usage to run symfc for fc3 with cutoff.
Total memory usage is memsize + memsize2. These are separated because
they behave differently with respect to cutoff distance.
batch_size is hardcoded to 100 because it is so in symfc.
"""
symfc_solver = SymfcFCSolver(supercell, symmetry, options={"cutoff": {3: cutoff}})
basis_size = symfc_solver.estimate_basis_size(orders=[3])[3]
memsize = basis_size**2 * 3 * 8 / 10**9
memsize2 = len(supercell) * 3 * batch_size * basis_size * 8 / 10**9
return memsize, memsize2
class FDFC3Solver:
"""Finite difference type force constants calculator.
@ -242,6 +111,193 @@ class FC3Solver(FCSolver):
def _set_traditional_solver(self, solver_class: Optional[type] = FDFC3Solver):
return super()._set_traditional_solver(solver_class=solver_class)
def _set_symfc_solver(self):
return super()._set_symfc_solver(order=3)
def _get_displacements_and_forces(self):
"""Return displacements and forces for fc3."""
return get_displacements_and_forces_fc3(self._dataset)
def extract_fc2_fc3_calculators(
fc_calculator: Optional[Union[str, dict]],
order: int,
) -> Optional[Union[str, dict]]:
"""Extract fc_calculator and fc_calculator_options for fc2 and fc3.
fc_calculator : str
FC calculator. "|" separates fc2 and fc3. First and last
parts separated correspond to fc2 and fc3 calculators, respectively.
order : int = 2 or 3
2 and 3 indicate fc2 and fc3, respectively.
"""
if isinstance(fc_calculator, dict) or fc_calculator is None:
return fc_calculator
elif isinstance(fc_calculator, str):
if "|" in fc_calculator:
_fc_calculator = fc_calculator.split("|")[order - 2]
if _fc_calculator == "":
return None
return _fc_calculator
else:
if fc_calculator.strip() == "":
return None
return fc_calculator
else:
raise RuntimeError("fc_calculator should be str, dict, or None.")
def update_cutoff_fc_calculator_options(
fc_calc_opts: Optional[Union[str, dict]],
cutoff_pair_distance: Optional[float],
) -> Optional[Union[str, dict]]:
"""Update fc_calculator_options with cutoff distances.
Parameters
----------
fc_calc_opts : str or dict
FC calculator options.
cutoff_pair_distance : float, optional
Cutoff distance for pair interaction.
"""
if cutoff_pair_distance is not None:
if not isinstance(fc_calc_opts, (str, dict)) and fc_calc_opts is not None:
raise RuntimeError("fc_calculator_options should be str, dict, or None.")
if isinstance(fc_calc_opts, dict) and "cutoff" not in fc_calc_opts:
fc_calc_opts["cutoff"] = float(cutoff_pair_distance)
elif isinstance(fc_calc_opts, str) and "cutoff" not in fc_calc_opts:
fc_calc_opts = f"{fc_calc_opts}, cutoff = {cutoff_pair_distance}"
elif fc_calc_opts is None:
fc_calc_opts = f"cutoff = {cutoff_pair_distance}"
return fc_calc_opts
def get_fc_calculator_params(
fc_calculator: Optional[str],
fc_calculator_options: Optional[str],
cutoff_pair_distance: Optional[float],
log_level: int = 0,
) -> tuple[Optional[str], Optional[str]]:
"""Compile fc_calculator and fc_calculator_options from input settings."""
_fc_calculator = None
fc_calculator_list = []
if fc_calculator is not None:
for fc_calculatr_str in fc_calculator.split("|"):
if fc_calculatr_str == "": # No external calculator
fc_calculator_list.append(fc_calculatr_str.lower())
elif fc_calculatr_str.lower() in fc_calculator_names:
fc_calculator_list.append(fc_calculatr_str.lower())
if fc_calculator_list:
_fc_calculator = "|".join(fc_calculator_list)
_fc_calculator_options = fc_calculator_options
if cutoff_pair_distance:
if fc_calculator_list and fc_calculator_list[-1] in ("alm", "symfc"):
if fc_calculator_list[-1] == "alm":
cutoff_str = f"-1 {cutoff_pair_distance}"
if fc_calculator_list[-1] == "symfc":
cutoff_str = f"{cutoff_pair_distance}"
_fc_calculator_options = _set_cutoff_in_fc_calculator_options(
_fc_calculator_options,
cutoff_str,
log_level,
)
return _fc_calculator, _fc_calculator_options
def determine_cutoff_pair_distance(
fc_calculator: Optional[str] = None,
fc_calculator_options: Optional[str] = None,
cutoff_pair_distance: Optional[float] = None,
symfc_memory_size: Optional[float] = None,
random_displacements: Optional[Union[int, str]] = None,
supercell: Optional[PhonopyAtoms] = None,
primitive: Optional[Primitive] = None,
symmetry: Optional[Symmetry] = None,
log_level: int = 0,
) -> float:
"""Determine cutoff pair distance for displacements."""
_cutoff_pair_distance, _symfc_memory_size = _get_cutoff_pair_distance(
fc_calculator,
fc_calculator_options,
cutoff_pair_distance,
symfc_memory_size,
)
if random_displacements is not None and random_displacements != "auto":
_symfc_memory_size = None
if _symfc_memory_size is not None:
if fc_calculator is None:
pass
elif fc_calculator != "symfc":
raise RuntimeError(
"Estimation of cutoff_pair_distance by memory size is only "
"available for symfc calculator."
)
symfc_options = {"memsize": {3: _symfc_memory_size}}
update_symfc_cutoff_by_memsize(
symfc_options, supercell, primitive, symmetry, verbose=log_level > 0
)
if symfc_options["cutoff"] is not None:
_cutoff_pair_distance = symfc_options["cutoff"][3]
return _cutoff_pair_distance
def _set_cutoff_in_fc_calculator_options(
fc_calculator_options: Optional[str],
cutoff_str: str,
log_level: int,
):
str_appended = f"cutoff={cutoff_str}"
calc_opts = fc_calculator_options
if calc_opts is None:
calc_opts = "|"
if "|" in calc_opts:
calc_opts_fc2, calc_opts_fc3 = [v.strip() for v in calc_opts.split("|")][:2]
else:
calc_opts_fc2 = calc_opts
calc_opts_fc3 = calc_opts
if calc_opts_fc3 == "":
calc_opts_fc3 += f"{str_appended}"
if log_level:
print(f'Set "{str_appended}" to fc_calculator_options for fc3.')
elif "cutoff" not in calc_opts_fc3:
calc_opts_fc3 += f", {str_appended}"
if log_level:
print(f'Appended "{str_appended}" to fc_calculator_options for fc3.')
return f"{calc_opts_fc2}|{calc_opts_fc3}"
def _get_cutoff_pair_distance(
fc_calculator: Optional[str],
fc_calculator_options: Optional[str],
cutoff_pair_distance: Optional[float],
symfc_memory_size: Optional[float] = None,
) -> Optional[float]:
"""Return cutoff_pair_distance from settings."""
_, _fc_calculator_options = get_fc_calculator_params(
fc_calculator,
fc_calculator_options,
cutoff_pair_distance,
)
symfc_options = parse_symfc_options(
extract_fc2_fc3_calculators(_fc_calculator_options, 3), 3
)
_cutoff_pair_distance = cutoff_pair_distance
cutoff = symfc_options.get("cutoff")
if cutoff is not None:
_cutoff_pair_distance = cutoff.get(3)
_symfc_memory_size = symfc_memory_size
memsize = symfc_options.get("memsize")
if memsize is not None:
_symfc_memory_size = memsize.get(3)
return _cutoff_pair_distance, _symfc_memory_size

164
phono3py/phonon/grid.py
View File

@ -39,12 +39,14 @@ from __future__ import annotations
import warnings
from collections.abc import Sequence
from numpy.typing import ArrayLike, NDArray
try:
from spglib import SpglibDataset
from spglib import SpglibDataset # type: ignore
except ImportError:
from types import SimpleNamespace as SpglibDataset
from types import SimpleNamespace
from typing import Optional, Union
import numpy as np
@ -132,6 +134,7 @@ class BZGrid:
D_diag : ndarray
P : ndarray
Q : ndarray
PS : ndarray
QDinv : ndarray
grid_matrix : ndarray
microzone_lattice : ndarray
@ -144,7 +147,7 @@ class BZGrid:
mesh: Union[int, float, Sequence, np.ndarray],
reciprocal_lattice=None,
lattice=None,
symmetry_dataset: Optional[Union[SpglibDataset]] = None,
symmetry_dataset: Optional[SpglibDataset] = None,
transformation_matrix: Optional[Union[Sequence, np.ndarray]] = None,
is_shift: Optional[Union[list, np.ndarray]] = None,
is_time_reversal: bool = True,
@ -202,17 +205,18 @@ class BZGrid:
self._is_shift = [v * 1 for v in is_shift]
self._is_time_reversal = is_time_reversal
self._store_dense_gp_map = store_dense_gp_map
self._addresses = None
self._gp_map = None
self._addresses: np.ndarray
self._gp_map: np.ndarray
self._grid_matrix = None
self._D_diag = np.ones(3, dtype="int64")
self._Q = np.eye(3, dtype="int64", order="C")
self._P = np.eye(3, dtype="int64", order="C")
self._QDinv = None
self._microzone_lattice = None
self._rotations = None
self._reciprocal_operations = None
self._gp_Gamma = None
self._QDinv: np.ndarray
self._microzone_lattice: np.ndarray
self._rotations: np.ndarray
self._reciprocal_operations: np.ndarray
self._rotations_cartesian: np.ndarray
self._gp_Gamma: int
if reciprocal_lattice is not None:
self._reciprocal_lattice = np.array(
@ -281,7 +285,7 @@ class BZGrid:
return np.array(np.dot(self.P, self._is_shift), dtype="int64")
@property
def grid_matrix(self):
def grid_matrix(self) -> Optional[np.ndarray]:
"""Grid generating matrix to be represented by SNF.
Grid generating matrix used for SNF.
@ -292,7 +296,7 @@ class BZGrid:
return self._grid_matrix
@property
def addresses(self):
def addresses(self) -> np.ndarray:
"""BZ-grid addresses.
Integer grid address of the points in Brillouin zone including
@ -315,7 +319,7 @@ class BZGrid:
return self._gp_map
@property
def gp_Gamma(self):
def gp_Gamma(self) -> int:
"""Return grid point index of Gamma-point."""
return self._gp_Gamma
@ -345,7 +349,7 @@ class BZGrid:
return self._grg2bzg
@property
def microzone_lattice(self):
def microzone_lattice(self) -> np.ndarray:
"""Basis vectors of microzone.
Basis vectors of microzone of GR-grid in column vectors.
@ -354,6 +358,16 @@ class BZGrid:
"""
return self._microzone_lattice
@property
def reciprocal_lattice(self) -> np.ndarray:
"""Reciprocal basis vectors of primitive cell.
Reciprocal basis vectors of primitive cell in column vectors.
shape=(3, 3), dtype='double', order='C'.
"""
return self._reciprocal_lattice
@property
def store_dense_gp_map(self):
"""Return gp_map type.
@ -364,7 +378,7 @@ class BZGrid:
return self._store_dense_gp_map
@property
def rotations(self):
def rotations(self) -> np.ndarray:
"""Return rotation matrices for grid points.
Rotation matrices for GR-grid addresses (g) defined as g'=Rg. This can
@ -376,12 +390,12 @@ class BZGrid:
return self._rotations
@property
def rotations_cartesian(self):
def rotations_cartesian(self) -> np.ndarray:
"""Return rotations in Cartesian coordinates."""
return self._rotations_cartesian
@property
def reciprocal_operations(self):
def reciprocal_operations(self) -> np.ndarray:
"""Return reciprocal rotations.
Reciprocal space rotation matrices in fractional coordinates defined as
@ -392,7 +406,7 @@ class BZGrid:
return self._reciprocal_operations
@property
def symmetry_dataset(self) -> SpglibDataset:
def symmetry_dataset(self) -> Optional[SpglibDataset]:
"""Return Symmetry.dataset."""
return self._symmetry_dataset
@ -539,7 +553,7 @@ class GridMatrix:
mesh: Union[int, float, Sequence, np.ndarray],
lattice: Union[Sequence, np.ndarray],
symmetry_dataset: Optional[SpglibDataset] = None,
transformation_matrix: Optional[Union[list, np.ndarray]] = None,
transformation_matrix: Optional[Union[Sequence, np.ndarray]] = None,
use_grg: bool = True,
force_SNF: bool = False,
SNF_coordinates: str = "reciprocal",
@ -591,7 +605,7 @@ class GridMatrix:
)
@property
def grid_matrix(self):
def grid_matrix(self) -> Optional[np.ndarray]:
"""Grid generating matrix to be represented by SNF.
Grid generating matrix used for SNF.
@ -602,7 +616,7 @@ class GridMatrix:
return self._grid_matrix
@property
def D_diag(self):
def D_diag(self) -> np.ndarray:
"""Diagonal elements of diagonal matrix after SNF: D=PAQ.
This corresponds to the mesh numbers in transformed reciprocal
@ -637,7 +651,7 @@ class GridMatrix:
mesh: Union[int, float, Sequence, np.ndarray],
use_grg: bool = False,
symmetry_dataset: Optional[SpglibDataset] = None,
transformation_matrix: Optional[Union[list, np.ndarray]] = None,
transformation_matrix: Optional[Union[Sequence, np.ndarray]] = None,
force_SNF=False,
coordinates="reciprocal",
) -> None:
@ -665,9 +679,8 @@ class GridMatrix:
have similar lengths.
"""
num_values = len(np.ravel(mesh))
if num_values == 1:
length = float(mesh)
try:
length = float(mesh) # type: ignore
if use_grg:
found_grg = self._run_grg(
symmetry_dataset,
@ -685,17 +698,20 @@ class GridMatrix:
length, self._lattice, rotations=symmetry_dataset.rotations
)
self._D_diag = np.array(mesh_numbers, dtype="int64")
if num_values == 9:
self._run_grg(
symmetry_dataset,
transformation_matrix,
None,
mesh,
force_SNF,
coordinates,
)
if num_values == 3:
self._D_diag = np.array(mesh, dtype="int64")
except (ValueError, TypeError):
num_values = len(np.ravel(mesh))
if num_values == 9:
self._run_grg(
symmetry_dataset,
transformation_matrix,
None,
mesh,
force_SNF,
coordinates,
)
if num_values == 3:
self._D_diag = np.array(mesh, dtype="int64")
def _run_grg(
self,
@ -730,7 +746,7 @@ class GridMatrix:
)
return False
def _get_mock_symmetry_dataset(self, transformation_matrix) -> dict:
def _get_mock_symmetry_dataset(self, transformation_matrix) -> SimpleNamespace:
"""Return mock symmetry_dataset containing transformation matrix.
Assuming self._lattice as standardized cell, and inverse of
@ -753,8 +769,6 @@ class GridMatrix:
)
raise RuntimeError(msg)
from types import SimpleNamespace
sym_dataset = SimpleNamespace(
**{
"rotations": np.eye(3, dtype="intc", order="C").reshape(1, 3, 3),
@ -768,7 +782,7 @@ class GridMatrix:
def _set_GRG_mesh(
self,
sym_dataset: SpglibDataset,
sym_dataset: Union[SpglibDataset, SimpleNamespace],
length: Optional[float] = None,
grid_matrix=None,
force_SNF=False,
@ -802,7 +816,10 @@ class GridMatrix:
self._grid_matrix = _grid_matrix
def _get_grid_matrix(
self, sym_dataset: dict, length: float, coordinates: str = "reciprocal"
self,
sym_dataset: Union[SpglibDataset, SimpleNamespace],
length: float,
coordinates: str = "reciprocal",
):
"""Return grid matrix.
@ -855,6 +872,22 @@ class GridMatrix:
return grid_matrix
def get_qpoints_from_bz_grid_points(
gps: Union[int, np.ndarray], bz_grid: BZGrid
) -> np.ndarray:
"""Return q-point(s) in reduced coordinates of grid point(s).
Parameters
----------
i_gps : int or ndarray
BZ-grid index (int) or indices (ndarray).
bz_grid : BZGrid
BZ-grid instance.
"""
return bz_grid.addresses[gps] @ bz_grid.QDinv.T
def get_grid_point_from_address_py(addresses, D_diag):
"""Return GR-grid point index from addresses.
@ -1041,12 +1074,12 @@ def _get_grid_points_by_bz_rotations_py(bz_gp, bz_grid: BZGrid, rotations):
)
+ num_grgp
).tolist()
gps.insert(0, gp)
gps.insert(0, gp) # type: ignore
indices = np.where((bz_grid.addresses[gps] == adrs).all(axis=1))[0]
if len(indices) == 0:
msg = "with_surface did not work properly."
raise RuntimeError(msg)
bzgps[i] = gps[indices[0]]
bzgps[i] = gps[indices[0]] # type: ignore
return bzgps
@ -1151,21 +1184,15 @@ def _relocate_BZ_grid_address(
else:
bz_map = np.zeros(np.prod(D_diag) * 9 + 1, dtype="int64")
# Mpr^-1 = Lr^-1 Lp
reclat_T = np.array(np.transpose(reciprocal_lattice), dtype="double", order="C")
reduced_basis = get_reduced_bases(reclat_T)
tmat_inv = np.dot(np.linalg.inv(reduced_basis.T), reclat_T.T)
tmat_inv_int = np.rint(tmat_inv).astype("int64")
assert (np.abs(tmat_inv - tmat_inv_int) < 1e-5).all()
reduced_basis, tmat_inv_int = get_reduced_bases_and_tmat_inv(reciprocal_lattice)
num_gp = recgrid.bz_grid_addresses(
bz_grid_addresses,
bz_map,
bzg2grg,
np.array(D_diag, dtype="int64"),
np.array(np.dot(tmat_inv_int, Q), dtype="int64", order="C"),
np.array(tmat_inv_int @ Q, dtype="int64", order="C"),
_PS,
np.array(reduced_basis.T, dtype="double", order="C"),
reduced_basis,
store_dense_gp_map * 1 + 1,
)
@ -1174,6 +1201,39 @@ def _relocate_BZ_grid_address(
return bz_grid_addresses, bz_map, bzg2grg
def get_reduced_bases_and_tmat_inv(
reciprocal_lattice: ArrayLike,
) -> tuple[NDArray, NDArray]:
"""Return reduced bases and inverse transformation matrix.
Parameters
----------
reciprocal_lattice : ArrayLike
Reciprocal lattice vectors in column vectors.
shape=(3, 3), dtype='double'
Returns
-------
reduced_basis : ndarray
Reduced basis vectors in column vectors.
shape=(3, 3), dtype='double', order='C'
tmat_inv_int : ndarray
Inverse transformation matrix in integer.
This is used to transform reciprocal lattice vectors to
conventional lattice vectors.
shape=(3, 3), dtype='int64'
"""
# Mpr^-1 = Lr^-1 Lp
reclat_T = np.array(np.transpose(reciprocal_lattice), dtype="double", order="C")
reduced_basis = get_reduced_bases(reclat_T)
assert reduced_basis is not None, "Reduced basis is not found."
tmat_inv = np.linalg.inv(reduced_basis.T) @ reclat_T.T
tmat_inv_int = np.rint(tmat_inv).astype("int64")
assert (np.abs(tmat_inv - tmat_inv_int) < 1e-5).all()
return np.array(reduced_basis.T, dtype="double", order="C"), tmat_inv_int
def _get_ir_grid_map(
D_diag: Union[np.ndarray, Sequence],
grg_rotations: Union[np.ndarray, Sequence],

331
phono3py/phonon3/interaction.py
View File

@ -35,11 +35,11 @@
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import warnings
from collections.abc import Sequence
from typing import Literal, Optional, Union
from typing import Literal
import numpy as np
from numpy.typing import NDArray
from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, get_dynamical_matrix
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import Primitive, compute_all_sg_permutations
@ -94,26 +94,26 @@ class Interaction:
primitive: Primitive,
bz_grid: BZGrid,
primitive_symmetry: Symmetry,
fc3: Optional[np.ndarray] = None,
band_indices: Optional[Union[np.ndarray, Sequence]] = None,
constant_averaged_interaction: Optional[float] = None,
frequency_factor_to_THz: Optional[float] = None,
frequency_scale_factor: Optional[float] = None,
unit_conversion: Optional[float] = None,
fc3: NDArray | None = None,
fc3_nonzero_indices: NDArray | None = None,
band_indices: NDArray | Sequence | None = None,
constant_averaged_interaction: float | None = None,
frequency_factor_to_THz: float | None = None,
frequency_scale_factor: float | None = None,
unit_conversion: float | None = None,
is_mesh_symmetry: bool = True,
symmetrize_fc3q: bool = False,
make_r0_average: bool = False,
cutoff_frequency: Optional[float] = None,
cutoff_frequency: float | None = None,
lapack_zheev_uplo: Literal["L", "U"] = "L",
openmp_per_triplets: Optional[bool] = None,
openmp_per_triplets: bool | None = None,
):
"""Init method."""
self._primitive = primitive
self._bz_grid = bz_grid
self._primitive_symmetry = primitive_symmetry
self._band_indices = None
self._set_band_indices(band_indices)
self._band_indices = self._get_band_indices(band_indices)
self._constant_averaged_interaction = constant_averaged_interaction
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().DefaultToTHz
@ -122,12 +122,12 @@ class Interaction:
self._frequency_scale_factor = frequency_scale_factor
if fc3 is not None:
self._set_fc3(fc3)
self._set_fc3(fc3, fc3_nonzero_indices=fc3_nonzero_indices)
# Unit to eV^2
if unit_conversion is None:
num_grid = np.prod(self.mesh_numbers)
self._unit_conversion = (
self._unit_conversion = float(
(get_physical_units().Hbar * get_physical_units().EV) ** 3
/ 36
/ 8
@ -147,7 +147,7 @@ class Interaction:
self._is_mesh_symmetry = is_mesh_symmetry
self._symmetrize_fc3q = symmetrize_fc3q
self._make_r0_average = make_r0_average
self._lapack_zheev_uplo = lapack_zheev_uplo
self._lapack_zheev_uplo: Literal["L", "U"] = lapack_zheev_uplo
self._openmp_per_triplets = openmp_per_triplets
self._symprec = self._primitive_symmetry.tolerance
@ -160,6 +160,7 @@ class Interaction:
self._g_zero = None
self._phonon_done = None
self._phonon_all_done = False
self._done_nac_at_gamma = False # Phonon at Gamma is calculated with NAC.
self._frequencies = None
self._eigenvectors = None
@ -181,13 +182,14 @@ class Interaction:
)
self._get_all_shortest()
def run(
self, lang: Literal["C", "Python"] = "C", g_zero: Optional[np.ndarray] = None
):
def run(self, lang: Literal["C", "Python"] = "C", g_zero: NDArray | None = None):
"""Run ph-ph interaction calculation."""
if (self._phonon_done == 0).any():
if self._phonon_all_done:
self.run_phonon_solver()
if self._triplets_at_q is None:
raise RuntimeError("Set grid point first by set_grid_point().")
num_band = len(self._primitive) * 3
num_triplets = len(self._triplets_at_q)
@ -207,7 +209,7 @@ class Interaction:
)
@property
def interaction_strength(self) -> np.ndarray:
def interaction_strength(self) -> NDArray | None:
"""Return ph-ph interaction strength.
Returns
@ -219,18 +221,8 @@ class Interaction:
"""
return self._interaction_strength
def get_interaction_strength(self):
"""Return ph-ph interaction strength."""
warnings.warn(
"Use attribute, Interaction.interaction_strength "
"instead of Interaction.get_interaction_strength().",
DeprecationWarning,
stacklevel=2,
)
return self.interaction_strength
@property
def mesh_numbers(self) -> np.ndarray:
def mesh_numbers(self) -> NDArray:
"""Return mesh numbers.
Returns
@ -241,73 +233,56 @@ class Interaction:
"""
return self._bz_grid.D_diag
def get_mesh_numbers(self):
"""Return mesh numbers."""
warnings.warn(
"Use attribute, Interaction.mesh_numbers "
"instead of Interaction.get_mesh_numbers().",
DeprecationWarning,
stacklevel=2,
)
return self.mesh_numbers
@property
def is_mesh_symmetry(self) -> bool:
"""Whether symmetry of grid is utilized or not."""
return self._is_mesh_symmetry
@property
def fc3(self) -> np.ndarray:
def fc3(self) -> NDArray:
"""Return fc3."""
return self._fc3
def get_fc3(self):
"""Return fc3."""
warnings.warn(
"Use attribute, Interaction.fc3 instead of Interaction.get_fc3().",
DeprecationWarning,
stacklevel=2,
)
return self.fc3
@property
def fc3_nonzero_indices(self) -> NDArray:
"""Return fc3_nonzero_indices."""
return self._fc3_nonzero_indices
@property
def dynamical_matrix(self) -> Optional[DynamicalMatrix]:
def dynamical_matrix(self) -> DynamicalMatrix | None:
"""Return DynamicalMatrix class instance."""
return self._dm
def get_dynamical_matrix(self):
"""Return DynamicalMatrix class instance."""
warnings.warn(
"Use attribute, Interaction.dynamical_matrix "
"instead of Interaction.get_dynamical_matrix().",
DeprecationWarning,
stacklevel=2,
)
return self.dynamical_matrix
@property
def primitive(self) -> Primitive:
"""Return Primitive class instance."""
return self._primitive
def get_primitive(self):
"""Return Primitive class instance."""
warnings.warn(
"Use attribute, Interaction.primitive "
"instead of Interaction.get_primitive().",
DeprecationWarning,
stacklevel=2,
)
return self.primitive
@property
def primitive_symmetry(self) -> Symmetry:
"""Return Symmetry class instance of primitive cell."""
return self._primitive_symmetry
@property
def phonon_all_done(self) -> bool:
"""Return whether phonon calculation is done at all grid points.
True value is set in run_phonon_solver(). Even when False, it is
possible that all phonons are already calculated, but it is safer to
think some of phonons are not calculated yet. To be sure, check
(self._phonon_done == 0).any().
"""
return self._phonon_all_done
def get_triplets_at_q(
self,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
) -> tuple[
NDArray | None,
NDArray | None,
NDArray | None,
NDArray | None,
]:
"""Return grid point triplets information.
triplets_at_q is in BZ-grid.
@ -329,7 +304,7 @@ class Interaction:
return self._bz_grid
@property
def band_indices(self) -> np.ndarray:
def band_indices(self) -> NDArray:
"""Return band indices.
Returns
@ -340,23 +315,13 @@ class Interaction:
"""
return self._band_indices
def get_band_indices(self):
"""Return band indices."""
warnings.warn(
"Use attribute, Interaction.band_indices "
"instead of Interaction.get_band_indices().",
DeprecationWarning,
stacklevel=2,
)
return self.band_indices
@property
def nac_params(self) -> dict:
def nac_params(self) -> dict | None:
"""Return NAC params."""
return self._nac_params
@property
def nac_q_direction(self) -> Optional[np.ndarray]:
def nac_q_direction(self) -> NDArray | None:
"""Return q-direction used for NAC at q->0.
Direction of q-vector watching from Gamma point used for
@ -376,28 +341,8 @@ class Interaction:
else:
self._nac_q_direction = np.array(nac_q_direction, copy=True, dtype="double")
def get_nac_q_direction(self):
"""Return q-direction used for NAC at q->0."""
warnings.warn(
"Use attribute, Interaction.nac_q_direction "
"instead of Interaction.get_nac_q_direction().",
DeprecationWarning,
stacklevel=2,
)
return self.nac_q_direction
def set_nac_q_direction(self, nac_q_direction=None):
"""Set NAC q-point direction valid at q->0."""
warnings.warn(
"Use attribute, Interaction.nac_q_direction "
"instead of Interaction.set_nac_q_direction().",
DeprecationWarning,
stacklevel=2,
)
self.nac_q_direction = nac_q_direction
@property
def zero_value_positions(self) -> Optional[np.ndarray]:
def zero_value_positions(self) -> NDArray | None:
"""Return zero ph-ph interaction elements information.
Returns
@ -407,17 +352,9 @@ class Interaction:
"""
return self._g_zero
def get_zero_value_positions(self):
"""Return zero ph-ph interaction elements information."""
warnings.warn(
"Use attribute, Interaction.zero_value_positions "
"instead of Interaction.get_zero_value_positions().",
DeprecationWarning,
stacklevel=2,
)
return self.zero_value_positions
def get_phonons(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
def get_phonons(
self,
) -> tuple[NDArray | None, NDArray | None, NDArray | None]:
"""Return phonons on grid.
Returns
@ -442,48 +379,18 @@ class Interaction:
"""Return phonon frequency conversion factor to THz."""
return self._frequency_factor_to_THz
def get_frequency_factor_to_THz(self):
"""Return phonon frequency conversion factor to THz."""
warnings.warn(
"Use attribute, Interaction.frequency_factor_to_THz ",
"instead of Interaction.get_frequency_factor_to_THz().",
DeprecationWarning,
stacklevel=2,
)
return self.frequency_factor_to_THz
@property
def lapack_zheev_uplo(self) -> Literal["L", "U"]:
"""Return U or L for lapack zheev solver."""
return self._lapack_zheev_uplo
def get_lapack_zheev_uplo(self):
"""Return U or L for lapack zheev solver."""
warnings.warn(
"Use attribute, Interaction.lapack_zheev_uplo "
"instead of Interaction.get_lapack_zheev_uplo().",
DeprecationWarning,
stacklevel=2,
)
return self.lapack_zheev_uplo
@property
def cutoff_frequency(self) -> float:
"""Return cutoff phonon frequency to judge imaginary phonon."""
return self._cutoff_frequency
def get_cutoff_frequency(self):
"""Return cutoff phonon frequency to judge imaginary phonon."""
warnings.warn(
"Use attribute, Interaction.cutoff_frequency "
"instead of Interaction.get_cutoff_frequency().",
DeprecationWarning,
stacklevel=2,
)
return self.cutoff_frequency
@property
def openmp_per_triplets(self) -> bool:
def openmp_per_triplets(self) -> bool | None:
"""Return whether OpenMP distribution over triplets or bands."""
return self._openmp_per_triplets
@ -505,7 +412,7 @@ class Interaction:
return self._make_r0_average
@property
def all_shortest(self) -> np.ndarray:
def all_shortest(self) -> NDArray:
"""Return boolean of make_r0_average.
This flag is used to activate averaging of fc3 transformation
@ -517,31 +424,27 @@ class Interaction:
return self._all_shortest
@property
def averaged_interaction(self) -> np.ndarray:
def averaged_interaction(self) -> NDArray:
"""Return sum over phonon triplets of interaction strength.
See Eq.(21) of PRB 91, 094306 (2015)
"""
if self._interaction_strength is None or self._weights_at_q is None:
raise RuntimeError(
"Interaction strength and weights at q are not set. "
"Run Interaction.run() first."
)
# v[triplet, band0, band, band]
v = self._interaction_strength
w = self._weights_at_q
v_sum = np.dot(w, v.sum(axis=2).sum(axis=2))
return v_sum / np.prod(v.shape[2:])
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,
stacklevel=2,
)
return self.averaged_interaction
def get_primitive_and_supercell_correspondence(
self,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
) -> tuple[NDArray, NDArray, NDArray, NDArray, NDArray]:
"""Return atomic pair information."""
return (self._svecs, self._multi, self._p2s, self._s2p, self._masses)
@ -550,31 +453,11 @@ class Interaction:
"""Return unit conversion factor."""
return self._unit_conversion
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,
stacklevel=2,
)
return self.unit_conversion_factor
@property
def constant_averaged_interaction(self) -> float:
def constant_averaged_interaction(self) -> float | None:
"""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,
stacklevel=2,
)
return self.constant_averaged_interaction
def set_interaction_strength(self, pp_strength, g_zero=None):
"""Set interaction strength."""
self._interaction_strength = pp_strength
@ -704,6 +587,14 @@ class Interaction:
"Input grid addresses are inconsistent. Setting phonons failed."
)
else:
if (
self._phonon_done is None
or self._frequencies is None
or self._eigenvectors is None
):
raise RuntimeError(
"Phonons are not initialized. Call init_dynamical_matrix() first."
)
self._phonon_done[:] = 1
self._frequencies[:] = frequencies
self._eigenvectors[:] = eigenvectors
@ -720,6 +611,7 @@ class Interaction:
self._run_phonon_solver_c(
np.arange(len(self._bz_grid.addresses), dtype="int64")
)
self._phonon_all_done = True
else:
self._run_phonon_solver_c(grid_points)
@ -737,6 +629,15 @@ class Interaction:
otherwise without NAC. Default is False.
"""
if (
self._phonon_done is None
or self._frequencies is None
or self._eigenvectors is None
):
raise RuntimeError(
"Phonons are not initialized. Call init_dynamical_matrix() first."
)
if not is_nac and self._frequencies_at_gamma is not None:
gp_Gamma = self._bz_grid.gp_Gamma
self._frequencies[gp_Gamma] = self._frequencies_at_gamma
@ -765,6 +666,11 @@ class Interaction:
self._phonon_done
"""
if self._phonon_done is None:
raise RuntimeError(
"Phonons are not initialized. Call init_dynamical_matrix() first."
)
self._phonon_done[:] = 0
ir_grid_points, _, _ = get_ir_grid_points(self._bz_grid)
ir_bz_grid_points = self._bz_grid.grg2bzg[ir_grid_points]
@ -775,8 +681,8 @@ class Interaction:
# perms.shape = (len(spg_ops), len(primitive)), dtype='intc'
perms = compute_all_sg_permutations(
self._primitive.scaled_positions,
self._bz_grid.symmetry_dataset.rotations,
self._bz_grid.symmetry_dataset.translations,
self._bz_grid.symmetry_dataset.rotations, # type: ignore
self._bz_grid.symmetry_dataset.translations, # type: ignore
np.array(self._primitive.cell.T, dtype="double", order="C"),
symprec=self._symprec,
)
@ -824,13 +730,13 @@ class Interaction:
"""
d2r_map = []
for r in self._bz_grid.symmetry_dataset.rotations:
for r in self._bz_grid.symmetry_dataset.rotations: # type: ignore
for i, rec_r in enumerate(self._bz_grid.reciprocal_operations):
if (rec_r.T == r).all():
d2r_map.append(i)
break
assert len(d2r_map) == len(self._bz_grid.symmetry_dataset.rotations)
assert len(d2r_map) == len(self._bz_grid.symmetry_dataset.rotations) # type: ignore
return d2r_map
@ -840,8 +746,12 @@ class Interaction:
e_j'(Rq) = R e_j(q) exp(-iRq.\tau)
"""
assert self._phonon_done is not None
assert self._frequencies is not None
assert self._eigenvectors is not None
Rq = np.dot(self._bz_grid.QDinv, self._bz_grid.addresses[bzgp])
tau = self._bz_grid.symmetry_dataset.translations[t_i]
tau = self._bz_grid.symmetry_dataset.translations[t_i] # type: ignore
phase_factor = np.exp(-2j * np.pi * np.dot(Rq, tau))
self._phonon_done[bzgp] = 1
self._frequencies[bzgp, :] = self._frequencies[orig_gp, :]
@ -863,6 +773,10 @@ class Interaction:
in this method.
"""
assert self._phonon_done is not None
assert self._frequencies is not None
assert self._eigenvectors is not None
r_inv = -np.eye(3, dtype="int64")
bz_grid_points_solved = []
for bzgp, done in enumerate(self._phonon_done):
@ -908,7 +822,7 @@ class Interaction:
self._interaction_strength = None
self._g_zero = None
def _set_fc3(self, fc3):
def _set_fc3(self, fc3: NDArray, fc3_nonzero_indices: NDArray | None = None):
if (
isinstance(fc3, np.ndarray)
and fc3.dtype == np.dtype("double")
@ -925,20 +839,42 @@ class Interaction:
fc3 * self._frequency_scale_factor**2, dtype="double", order="C"
)
def _set_band_indices(self, band_indices):
if fc3_nonzero_indices is None:
self._fc3_nonzero_indices = np.ones(
self._fc3.shape[:3], dtype="byte", order="C"
)
elif (
isinstance(fc3_nonzero_indices, np.ndarray)
and fc3_nonzero_indices.dtype == np.dtype("byte")
and fc3_nonzero_indices.flags.aligned
and fc3_nonzero_indices.flags.owndata
and fc3_nonzero_indices.flags.c_contiguous
):
self._fc3_nonzero_indices = fc3_nonzero_indices
else:
self._fc3_nonzero_indices = np.array(
fc3_nonzero_indices, dtype="byte", order="C"
)
def _get_band_indices(self, band_indices) -> NDArray:
num_band = len(self._primitive) * 3
if band_indices is None:
self._band_indices = np.arange(num_band, dtype="int64")
return np.arange(num_band, dtype="int64")
else:
self._band_indices = np.array(band_indices, dtype="int64")
return np.array(band_indices, dtype="int64")
def _run_c(self, g_zero):
import phono3py._phono3py as phono3c
assert self._interaction_strength is not None
assert self._triplets_at_q is not None
num_band = len(self._primitive) * 3
if g_zero is None or self._symmetrize_fc3q:
_g_zero = np.zeros(
self._interaction_strength.shape, dtype="byte", order="C"
self._interaction_strength.shape,
dtype="byte",
order="C",
)
else:
_g_zero = g_zero
@ -963,6 +899,7 @@ class Interaction:
self._bz_grid.D_diag,
self._bz_grid.Q,
self._fc3,
self._fc3_nonzero_indices,
self._svecs,
self._multi,
self._masses,
@ -993,6 +930,11 @@ class Interaction:
)
def _run_py(self):
assert self._interaction_strength is not None
assert self._triplets_at_q is not None
assert self._frequencies is not None
assert self._eigenvectors is not None
r2r = RealToReciprocal(
self._fc3, self._primitive, self.mesh_numbers, symprec=self._symprec
)
@ -1011,8 +953,10 @@ class Interaction:
for gp in grid_triplet:
self._run_phonon_solver_py(gp)
r2n.run(fc3_reciprocal, grid_triplet)
fc3_normal = r2n.get_reciprocal_to_normal()
assert fc3_normal is not None
self._interaction_strength[i] = (
np.abs(r2n.get_reciprocal_to_normal()) ** 2 * self._unit_conversion
np.abs(fc3_normal) ** 2 * self._unit_conversion
)
def _run_phonon_solver_py(self, grid_point):
@ -1038,6 +982,7 @@ class Interaction:
num_band = len(self._primitive) * 3
num_grid = len(self._bz_grid.addresses)
self._phonon_done = np.zeros(num_grid, dtype="byte")
self._phonon_all_done = False
self._frequencies = np.zeros((num_grid, num_band), dtype="double", order="C")
complex_dtype = "c%d" % (np.dtype("double").itemsize * 2)
self._eigenvectors = np.zeros(

21
phono3py/phonon3/reciprocal_to_normal.py
View File

@ -34,7 +34,11 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from collections.abc import Sequence
from typing import Optional, Union
import numpy as np
from phonopy.structure.atoms import PhonopyAtoms
class ReciprocalToNormal:
@ -47,7 +51,12 @@ class ReciprocalToNormal:
"""
def __init__(
self, primitive, frequencies, eigenvectors, band_indices, cutoff_frequency=0
self,
primitive: PhonopyAtoms,
frequencies: np.ndarray,
eigenvectors: np.ndarray,
band_indices: Union[Sequence[int], np.ndarray],
cutoff_frequency: float = 0,
):
"""Init method."""
self._primitive = primitive
@ -56,10 +65,8 @@ class ReciprocalToNormal:
self._band_indices = band_indices
self._cutoff_frequency = cutoff_frequency
self._masses = self._primitive.masses
self._fc3_normal = None
self._fc3_reciprocal = None
self._fc3_normal: np.ndarray
self._fc3_reciprocal: np.ndarray
def run(self, fc3_reciprocal, grid_triplet):
"""Calculate fc3 in phonon coordinates."""
@ -71,7 +78,7 @@ class ReciprocalToNormal:
)
self._reciprocal_to_normal(grid_triplet)
def get_reciprocal_to_normal(self):
def get_reciprocal_to_normal(self) -> Optional[np.ndarray]:
"""Return fc3 in phonon coordinates."""
return self._fc3_normal
@ -102,7 +109,7 @@ class ReciprocalToNormal:
* e3[k * 3 + n, b3]
* self._fc3_reciprocal[i, j, k, ll, m, n]
)
mass_sqrt = np.sqrt(np.prod(self._masses[[i, j, k]]))
mass_sqrt = np.sqrt(np.prod(self._primitive.masses[[i, j, k]]))
sum_fc3 += sum_fc3_cart / mass_sqrt
return sum_fc3

14
phono3py/phonon3/triplets.py
View File

@ -42,7 +42,11 @@ from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phono3py.other.tetrahedron_method import get_tetrahedra_relative_grid_address
from phono3py.phonon.func import gaussian
from phono3py.phonon.grid import BZGrid, get_grid_point_from_address_py
from phono3py.phonon.grid import (
BZGrid,
get_grid_point_from_address_py,
get_reduced_bases_and_tmat_inv,
)
if TYPE_CHECKING:
from phono3py.phonon3.interaction import Interaction
@ -327,14 +331,18 @@ def _get_BZ_triplets_at_q(bz_grid_index, bz_grid: BZGrid, map_triplets):
weights[g] += 1
ir_weights = np.extract(weights > 0, weights)
triplets = -np.ones((len(ir_weights), 3), dtype="int64", order="C")
reduced_basis, tmat_inv_int = get_reduced_bases_and_tmat_inv(
bz_grid.reciprocal_lattice
)
num_ir_ret = phono3c.BZ_triplets_at_q(
triplets,
bz_grid_index,
bz_grid.addresses,
bz_grid.gp_map,
map_triplets,
np.array(bz_grid.D_diag, dtype="int64"),
bz_grid.Q,
bz_grid.D_diag,
np.array(tmat_inv_int @ bz_grid.Q, dtype="int64", order="C"),
reduced_basis,
bz_grid.store_dense_gp_map * 1 + 1,
)

2
phono3py/version.py
View File

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

2
pyproject.toml
View File

@ -16,7 +16,7 @@ dependencies = [
"matplotlib",
"h5py",
"spglib",
"phonopy>=2.38.2,<2.39",
"phonopy>=2.39,<2.40",
]
license = "BSD-3-Clause"
license-files = ["LICENSE"]

0
test/cui/test_phono3py_load.py → test/api/test_phono3py_load.py
View File

12
test/conductivity/test_kappa_RTA.py
View File

@ -276,6 +276,18 @@ def test_kappa_RTA_nacl(nacl_pbe: Phono3py):
np.testing.assert_allclose(ref_kappa_RTA, kappa, atol=0.5)
def test_kappa_RTA_mgo(
mgo_222rd_444rd_symfc: Phono3py, mgo_222rd_444rd_symfc_compact_fc: Phono3py
):
"""Test RTA by MgO cutoff 4."""
for ph3 in (mgo_222rd_444rd_symfc, mgo_222rd_444rd_symfc_compact_fc):
ref_kappa_RTA = [63.75, 63.75, 63.75, 0, 0, 0]
kappa = _get_kappa(ph3, [11, 11, 11]).ravel()
np.testing.assert_allclose(ref_kappa_RTA, kappa, atol=0.5)
kappa = _get_kappa(ph3, [11, 11, 11], is_full_pp=True).ravel()
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:

85
test/conftest.py
View File

@ -322,6 +322,49 @@ def si_pbesol_111_222_alm_cutoff(request) -> Phono3py:
)
@pytest.fixture(scope="session")
def si_pbesol_111_222_symfc_cutoff() -> Phono3py:
"""Return Phono3py instance of Si 1x1x1.
* with symmetry
* full fc
* use symfc if available on test side
* cutoff=3
"""
pytest.importorskip("symfc")
yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
return phono3py.load(
yaml_filename,
fc_calculator="symfc",
fc_calculator_options="cutoff = 3",
log_level=1,
)
@pytest.fixture(scope="session")
def si_pbesol_111_222_symfc_cutoff_compact_fc() -> Phono3py:
"""Return Phono3py instance of Si 1x1x1.
* with symmetry
* compact fc
* use symfc if available on test side
* cutoff=3
"""
pytest.importorskip("symfc")
yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
return phono3py.load(
yaml_filename,
fc_calculator="symfc",
fc_calculator_options="cutoff = 3",
is_compact_fc=True,
log_level=1,
)
@pytest.fixture(scope="session")
def si_pbesol_111_222_alm_cutoff_fc2(request) -> Phono3py:
"""Return Phono3py instance of Si 1x1x1.
@ -523,6 +566,48 @@ def mgo_222rd_444rd() -> Phono3py:
return phono3py.load(yaml_filename, produce_fc=False, log_level=1)
@pytest.fixture(scope="session")
def mgo_222rd_444rd_symfc() -> Phono3py:
"""Return Phono3py instance of MgO-2x2x2-4x4x4 RD-RD.
* with symmetry
* full fc
* use symfc if available on test side
"""
pytest.importorskip("symfc")
yaml_filename = cwd / "phono3py_params_MgO-222rd-444rd.yaml.xz"
return phono3py.load(
yaml_filename,
is_compact_fc=False,
fc_calculator="symfc",
fc_calculator_options="|cutoff = 4",
log_level=1,
)
@pytest.fixture(scope="session")
def mgo_222rd_444rd_symfc_compact_fc() -> Phono3py:
"""Return Phono3py instance of MgO-2x2x2-4x4x4 RD-RD.
* with symmetry
* full fc
* use symfc if available on test side
"""
pytest.importorskip("symfc")
yaml_filename = cwd / "phono3py_params_MgO-222rd-444rd.yaml.xz"
return phono3py.load(
yaml_filename,
is_compact_fc=True,
fc_calculator="symfc",
fc_calculator_options="|cutoff = 4",
log_level=1,
)
@pytest.fixture(scope="session")
def ph_nacl() -> Phonopy:
"""Return Phonopy class instance of NaCl 2x2x2."""

262
test/cui/test_phono3py_load_script.py
View File

@ -12,6 +12,7 @@ import h5py
import numpy as np
import pytest
import phono3py
from phono3py.cui.phono3py_script import main
cwd = pathlib.Path(__file__).parent
@ -22,24 +23,27 @@ cwd_called = pathlib.Path.cwd()
class MockArgs:
"""Mock args of ArgumentParser."""
filename: Optional[Sequence[os.PathLike]] = None
cell_filename: Optional[str] = None
conf_filename: Optional[os.PathLike] = None
fc_calculator: Optional[str] = None
fc_calculator_options: Optional[str] = None
fc_symmetry: bool = True
filename: Optional[Sequence[os.PathLike]] = None
force_sets_mode: bool = False
force_sets_to_forces_fc2_mode: bool = False
input_filename = None
input_output_filename = None
log_level: Optional[int] = None
output_yaml_filename: Optional[os.PathLike] = None
show_num_triplets: bool = False
write_grid_points: bool = False
fc_symmetry: bool = True
cell_filename: Optional[str] = None
is_bterta: Optional[bool] = None
mesh_numbers: Optional[Sequence] = None
mlp_params: Optional[str] = None
output_filename = None
output_yaml_filename: Optional[os.PathLike] = None
random_displacements: Optional[Union[int, str]] = None
show_num_triplets: bool = False
temperatures: Optional[Sequence] = None
use_pypolymlp: bool = False
input_filename = None
output_filename = None
input_output_filename = None
write_grid_points: bool = False
def __iter__(self):
"""Make self iterable to support in."""
@ -54,7 +58,61 @@ def test_phono3py_load():
"""Test phono3py-load script."""
# Check sys.exit(0)
argparse_control = _get_phono3py_load_args(
cwd / ".." / "phono3py_params_Si-111-222.yaml"
cwd / ".." / "phono3py_params_Si-111-222.yaml",
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
assert excinfo.value.code == 0
argparse_control = _get_phono3py_load_args(
cwd_called / "phono3py.yaml",
is_bterta=True,
temperatures=[
"300",
],
mesh_numbers=["5", "5", "5"],
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
assert excinfo.value.code == 0
# Clean files created by phono3py-load script.
for created_filename in (
"phono3py.yaml",
"fc2.hdf5",
"fc3.hdf5",
"kappa-m555.hdf5",
):
file_path = cwd_called / created_filename
if file_path.exists():
file_path.unlink()
@pytest.mark.parametrize(
"load_phono3py_yaml,fc_calculator,fc_calculator_options",
[
(True, None, None),
(True, "symfc", None),
(True, "symfc", "|cutoff=4.0"),
(False, "symfc", "|cutoff=4.0"),
],
)
def test_phono3py_load_with_typeII_dataset(
fc_calculator: str | None,
fc_calculator_options: str | None,
load_phono3py_yaml: bool,
):
"""Test phono3py-load script with typeII dataset.
When None, fallback to symfc.
"""
pytest.importorskip("symfc")
argparse_control = _get_phono3py_load_args(
cwd / ".." / "phono3py_params-Si111-rd.yaml.xz",
load_phono3py_yaml=load_phono3py_yaml,
fc_calculator=fc_calculator,
fc_calculator_options=fc_calculator_options,
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
@ -64,28 +122,14 @@ def test_phono3py_load():
for created_filename in ("phono3py.yaml", "fc2.hdf5", "fc3.hdf5"):
file_path = pathlib.Path(cwd_called / created_filename)
if file_path.exists():
file_path.unlink()
@pytest.mark.parametrize("fc_calculator,exit_code", [(None, 0), ("symfc", 0)])
def test_phono3py_load_with_typeII_dataset(fc_calculator, exit_code):
"""Test phono3py-load script with typeII dataset.
When None, fallback to symfc.
"""
pytest.importorskip("symfc")
argparse_control = _get_phono3py_load_args(
cwd / ".." / "phono3py_params-Si111-rd.yaml.xz", fc_calculator=fc_calculator
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
assert excinfo.value.code == exit_code
# Clean files created by phono3py-load script.
for created_filename in ("phono3py.yaml", "fc2.hdf5", "fc3.hdf5"):
file_path = pathlib.Path(cwd_called / created_filename)
if file_path.exists():
if created_filename == "fc3.hdf5":
with h5py.File(file_path, "r") as f:
if fc_calculator_options is None:
assert "fc3_nonzero_indices" not in f
else:
assert "fc3_nonzero_indices" in f
assert "fc3_cutoff" in f
assert f["fc3_cutoff"][()] == pytest.approx(4.0)
file_path.unlink()
@ -129,25 +173,36 @@ def test_phono3py_load_with_pypolymlp_si():
pytest.importorskip("pypolymlp", minversion="0.9.2")
pytest.importorskip("symfc")
# Create fc2.hdf5
argparse_control = _get_phono3py_load_args(
cwd / ".." / "phono3py_params_Si-111-222-rd.yaml.xz",
fc_calculator="symfc",
use_pypolymlp=True,
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
assert excinfo.value.code == 0
# phono3py.yaml and fc2.hd5 are used in the next run. So they are not deleted.
for created_filename in ("fc3.hdf5", "phono3py_mlp_eval_dataset.yaml"):
for created_filename in ("phono3py.yaml", "fc2.hdf5", "fc3.hdf5"):
file_path = pathlib.Path(cwd_called / created_filename)
if file_path.exists():
file_path.unlink()
assert file_path.exists()
pathlib.Path(cwd_called / "fc3.hdf5").unlink()
# Create MLP (polymlp.yaml)
argparse_control = _get_phono3py_load_args(
cwd / ".." / "phono3py_params_Si-111-222-rd.yaml.xz",
use_pypolymlp=True,
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
assert excinfo.value.code == 0
for created_filename in ("phono3py.yaml", "polymlp.yaml"):
file_path = pathlib.Path(cwd_called / created_filename)
assert file_path.exists()
# Create phono3py_mlp_eval_dataset.yaml
argparse_control = _get_phono3py_load_args(
cwd_called / "phono3py.yaml",
fc_calculator="symfc",
random_displacements="auto",
use_pypolymlp=True,
)
@ -155,6 +210,9 @@ def test_phono3py_load_with_pypolymlp_si():
main(**argparse_control)
assert excinfo.value.code == 0
ph3 = phono3py.load(cwd_called / "phono3py_mlp_eval_dataset.yaml")
assert len(ph3.displacements) == 4
for created_filename in (
"phono3py.yaml",
"fc2.hdf5",
@ -163,17 +221,123 @@ def test_phono3py_load_with_pypolymlp_si():
"phono3py_mlp_eval_dataset.yaml",
):
file_path = pathlib.Path(cwd_called / created_filename)
if file_path.exists():
file_path.unlink()
assert file_path.exists()
file_path.unlink()
def test_phono3py_load_with_pypolymlp_nacl():
"""Test phono3py-load script with pypolymlp using NaCl.
First run generates polymlp.yaml.
Second run uses polymlp.yaml.
"""
pytest.importorskip("pypolymlp", minversion="0.9.2")
pytest.importorskip("symfc")
# Stage1 (preparation)
argparse_control = _get_phono3py_load_args(
cwd / ".." / "phono3py_params_MgO-222rd-444rd.yaml.xz",
mlp_params="cutoff=4.0,gtinv_maxl=4 4,max_p=1,gtinv_order=2",
fc_calculator="symfc",
random_displacements="auto",
use_pypolymlp=True,
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
assert excinfo.value.code == 0
ph3 = phono3py.load(cwd_called / "phono3py_mlp_eval_dataset.yaml")
assert len(ph3.displacements) == 16
for created_filename in (
"phono3py.yaml",
"fc2.hdf5",
"fc3.hdf5",
"polymlp.yaml",
"phono3py_mlp_eval_dataset.yaml",
):
file_path = pathlib.Path(cwd_called / created_filename)
assert file_path.exists()
for created_filename in (
"fc3.hdf5",
"phono3py_mlp_eval_dataset.yaml",
):
file_path = pathlib.Path(cwd_called / created_filename)
assert file_path.exists()
file_path.unlink()
# Stage2 (cutoff test)
argparse_control = _get_phono3py_load_args(
cwd_called / "phono3py.yaml",
fc_calculator="symfc",
fc_calculator_options="|cutoff=4.0",
random_displacements="auto",
use_pypolymlp=True,
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
assert excinfo.value.code == 0
ph3 = phono3py.load(cwd_called / "phono3py_mlp_eval_dataset.yaml")
assert len(ph3.displacements) == 4
for created_filename in (
"phono3py.yaml",
"fc2.hdf5",
"fc3.hdf5",
"polymlp.yaml",
"phono3py_mlp_eval_dataset.yaml",
):
file_path = pathlib.Path(cwd_called / created_filename)
assert file_path.exists()
for created_filename in (
"fc3.hdf5",
"phono3py_mlp_eval_dataset.yaml",
):
file_path = pathlib.Path(cwd_called / created_filename)
assert file_path.exists()
file_path.unlink()
# Stage3 (memsize test)
argparse_control = _get_phono3py_load_args(
cwd_called / "phono3py.yaml",
fc_calculator="symfc",
fc_calculator_options="|memsize=0.05",
random_displacements="auto",
use_pypolymlp=True,
)
with pytest.raises(SystemExit) as excinfo:
main(**argparse_control)
assert excinfo.value.code == 0
ph3 = phono3py.load(cwd_called / "phono3py_mlp_eval_dataset.yaml")
assert len(ph3.displacements) == 8
for created_filename in (
"phono3py.yaml",
"fc2.hdf5",
"fc3.hdf5",
"polymlp.yaml",
"phono3py_mlp_eval_dataset.yaml",
):
file_path = pathlib.Path(cwd_called / created_filename)
assert file_path.exists()
file_path.unlink()
def _get_phono3py_load_args(
phono3py_yaml_filepath: Union[str, pathlib.Path],
fc_calculator: Optional[str] = None,
fc_calculator_options: Optional[str] = None,
load_phono3py_yaml: bool = True,
is_bterta: bool = False,
temperatures: Optional[Sequence] = None,
mesh_numbers: Optional[Sequence] = None,
mlp_params: Optional[str] = None,
random_displacements: Optional[Union[int, str]] = None,
temperatures: Optional[Sequence] = None,
use_pypolymlp: bool = False,
):
# Mock of ArgumentParser.args.
@ -181,21 +345,27 @@ def _get_phono3py_load_args(
mockargs = MockArgs(
filename=[phono3py_yaml_filepath],
fc_calculator=fc_calculator,
fc_calculator_options=fc_calculator_options,
is_bterta=is_bterta,
temperatures=temperatures,
mesh_numbers=mesh_numbers,
use_pypolymlp=use_pypolymlp,
log_level=1,
mesh_numbers=mesh_numbers,
mlp_params=mlp_params,
random_displacements=random_displacements,
temperatures=temperatures,
use_pypolymlp=use_pypolymlp,
)
else:
mockargs = MockArgs(
filename=[],
fc_calculator=fc_calculator,
fc_calculator_options=fc_calculator_options,
log_level=1,
cell_filename=phono3py_yaml_filepath,
is_bterta=is_bterta,
temperatures=temperatures,
mesh_numbers=mesh_numbers,
mlp_params=mlp_params,
random_displacements=random_displacements,
temperatures=temperatures,
use_pypolymlp=use_pypolymlp,
)

69
test/interface/test_fc_calculator.py Normal file
View File

@ -0,0 +1,69 @@
"""Tests of functions in fc_calculator."""
import pytest
from phonopy.structure.atoms import PhonopyAtoms
from phono3py.api_phono3py import Phono3py
from phono3py.interface.fc_calculator import determine_cutoff_pair_distance
def test_determine_cutoff_pair_distance() -> None:
"""Test determine_cutoff_pair_distance."""
cutoff = determine_cutoff_pair_distance(fc_calculator_options="|cutoff=4")
assert cutoff == pytest.approx(4.0)
cutoff = determine_cutoff_pair_distance(cutoff_pair_distance=5.0)
assert cutoff == pytest.approx(5.0)
cutoff = determine_cutoff_pair_distance(
fc_calculator_options="|cutoff=4", cutoff_pair_distance=5.0
)
assert cutoff == pytest.approx(4.0)
def test_determine_cutoff_pair_distance_with_memsize(aln_cell: PhonopyAtoms) -> None:
"""Test determine_cutoff_pair_distance estimated by memsize."""
pytest.importorskip("symfc")
ph3 = Phono3py(aln_cell, supercell_matrix=[3, 3, 2])
cutoff = determine_cutoff_pair_distance(
fc_calculator="symfc",
fc_calculator_options="|memsize=0.1",
supercell=ph3.supercell,
primitive=ph3.primitive,
symmetry=ph3.symmetry,
log_level=1,
)
assert cutoff == pytest.approx(3.2)
cutoff = determine_cutoff_pair_distance(
fc_calculator="symfc",
symfc_memory_size=0.2,
supercell=ph3.supercell,
primitive=ph3.primitive,
symmetry=ph3.symmetry,
log_level=1,
)
assert cutoff == pytest.approx(3.7)
cutoff = determine_cutoff_pair_distance(
fc_calculator="symfc",
fc_calculator_options="|memsize=0.1",
symfc_memory_size=0.2,
supercell=ph3.supercell,
primitive=ph3.primitive,
symmetry=ph3.symmetry,
log_level=1,
)
assert cutoff == pytest.approx(3.2)
with pytest.raises(RuntimeError):
cutoff = determine_cutoff_pair_distance(
fc_calculator="alm",
fc_calculator_options="|memsize=0.1",
symfc_memory_size=0.2,
supercell=ph3.supercell,
primitive=ph3.primitive,
symmetry=ph3.symmetry,
log_level=1,
)

30
test/phonon3/test_fc3.py
View File

@ -218,10 +218,40 @@ def test_phonon_smat_fd_symfc(si_pbesol_111_222_fd_symfc: Phono3py):
def test_phonon_smat_alm_cutoff(si_pbesol_111_222_alm_cutoff: Phono3py):
"""Test phonon smat and alm with Si PBEsol 1x1x1-2x2x2 cutoff."""
ph = si_pbesol_111_222_alm_cutoff
assert ph.fc3 is not None
assert ph.fc2 is not None
np.testing.assert_allclose(ph.fc3[0, 1, 7], 0, atol=1e-6, rtol=0)
np.testing.assert_allclose(ph.fc2[0, 33], 0, atol=1e-6, rtol=0)
def test_phonon_smat_symfc_cutoff(si_pbesol_111_222_symfc_cutoff: Phono3py):
"""Test phonon smat and symfc with Si PBEsol 1x1x1-2x2x2 cutoff."""
ph = si_pbesol_111_222_symfc_cutoff
assert ph.fc3 is not None
assert ph.fc2 is not None
assert ph.fc3_nonzero_indices is not None
np.testing.assert_allclose(ph.fc3[0, 1, 7], 0, atol=1e-6, rtol=0)
np.testing.assert_allclose(ph.fc2[0, 33], 0, atol=1e-6, rtol=0)
assert ph.fc3.shape == (8, 8, 8, 3, 3, 3)
assert ph.fc3_nonzero_indices.shape == (8, 8, 8)
assert ph.fc3_nonzero_indices.sum() == 104
def test_phonon_smat_symfc_cutoff_compact_fc(
si_pbesol_111_222_symfc_cutoff_compact_fc: Phono3py,
):
"""Test phonon smat and symfc with Si PBEsol 1x1x1-2x2x2 cutoff."""
ph = si_pbesol_111_222_symfc_cutoff_compact_fc
assert ph.fc3 is not None
assert ph.fc2 is not None
assert ph.fc3_nonzero_indices is not None
np.testing.assert_allclose(ph.fc3[0, 1, 7], 0, atol=1e-6, rtol=0)
np.testing.assert_allclose(ph.fc2[0, 33], 0, atol=1e-6, rtol=0)
assert ph.fc3.shape == (2, 8, 8, 3, 3, 3)
assert ph.fc3_nonzero_indices.shape == (2, 8, 8)
assert ph.fc3_nonzero_indices.sum() == 26
def test_phonon_smat_alm_cutoff_fc2(si_pbesol_111_222_alm_cutoff_fc2: Phono3py):
"""Test phonon smat and alm with Si PBEsol 1x1x1-2x2x2 cutoff fc2."""
ph = si_pbesol_111_222_alm_cutoff_fc2

2
test/phonon3/test_interaction.py
View File

@ -301,6 +301,7 @@ def test_get_all_shortest(aln_lda: Phono3py):
"""Test Interaction._get_all_shortest."""
ph3 = aln_lda
ph3.mesh_numbers = 30
assert ph3.grid is not None
itr = Interaction(
ph3.primitive,
ph3.grid,
@ -339,6 +340,7 @@ def _get_irt(
make_r0_average: bool = False,
):
ph3.mesh_numbers = mesh
assert ph3.grid is not None
itr = Interaction(
ph3.primitive,
ph3.grid,

97
test/phonon3/test_triplets.py
View File

@ -911,7 +911,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[1, 36, 60],
[1, 38, 59],
[1, 41, 56],
],
], # 0
[
[1, 0, 4],
[1, 1, 3],
@ -937,7 +937,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[1, 36, 60],
[1, 38, 59],
[1, 41, 56],
],
], # 1
[
[1, 0, 4],
[1, 1, 3],
@ -969,7 +969,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[1, 39, 57],
[1, 41, 56],
[1, 42, 55],
],
], # 2
[
[1, 0, 4],
[1, 1, 3],
@ -1016,7 +1016,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[1, 70, 26],
[1, 72, 25],
[1, 73, 24],
],
], # 3
[
[1, 0, 4],
[1, 1, 3],
@ -1150,11 +1150,11 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[8, 2, 87],
[8, 3, 86],
[9, 4, 92],
[8, 5, 84],
[8, 6, 82],
[8, 8, 81],
[8, 10, 80],
[8, 11, 85],
[9, 5, 84],
[8, 7, 82],
[8, 9, 81],
[9, 10, 80],
[9, 11, 85],
[8, 12, 78],
[8, 13, 76],
[8, 14, 75],
@ -1164,23 +1164,23 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[9, 22, 67],
[8, 24, 65],
[8, 27, 62],
[8, 29, 66],
[9, 30, 66],
[8, 31, 58],
[8, 32, 57],
[8, 40, 50],
[8, 48, 48],
],
], # 0
[
[8, 0, 89],
[8, 1, 88],
[8, 2, 87],
[8, 3, 86],
[9, 4, 92],
[8, 5, 84],
[8, 6, 82],
[8, 8, 81],
[8, 10, 80],
[8, 11, 85],
[9, 5, 84],
[8, 7, 82],
[8, 9, 81],
[9, 10, 80],
[9, 11, 85],
[8, 12, 78],
[8, 13, 76],
[8, 14, 75],
@ -1190,23 +1190,23 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[9, 22, 67],
[8, 24, 65],
[8, 27, 62],
[8, 29, 66],
[9, 30, 66],
[8, 31, 58],
[8, 32, 57],
[8, 40, 50],
[8, 48, 48],
],
], # 1
[
[8, 0, 89],
[8, 1, 88],
[8, 2, 87],
[8, 3, 86],
[9, 4, 92],
[8, 5, 84],
[8, 6, 82],
[8, 8, 81],
[8, 10, 80],
[8, 11, 85],
[9, 5, 84],
[8, 7, 82],
[8, 9, 81],
[9, 10, 80],
[9, 11, 85],
[8, 12, 78],
[8, 13, 76],
[8, 14, 75],
@ -1219,7 +1219,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[8, 24, 65],
[8, 25, 64],
[8, 27, 62],
[8, 29, 66],
[9, 30, 66],
[8, 31, 58],
[8, 32, 57],
[8, 33, 56],
@ -1227,7 +1227,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[8, 40, 50],
[9, 41, 49],
[9, 42, 54],
[8, 43, 46],
[8, 43, 47],
[8, 44, 45],
[8, 48, 48],
[8, 50, 40],
@ -1240,20 +1240,20 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[8, 71, 19],
[9, 72, 18],
[8, 79, 17],
[8, 81, 8],
[8, 81, 9],
[8, 89, 0],
],
], # 2
[
[8, 0, 89],
[8, 1, 88],
[8, 2, 87],
[8, 3, 86],
[9, 4, 92],
[8, 5, 84],
[8, 6, 82],
[8, 8, 81],
[8, 10, 80],
[8, 11, 85],
[9, 5, 84],
[8, 7, 82],
[8, 9, 81],
[9, 10, 80],
[9, 11, 85],
[8, 12, 78],
[8, 13, 76],
[8, 14, 75],
@ -1266,7 +1266,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[8, 24, 65],
[8, 25, 64],
[8, 27, 62],
[8, 29, 66],
[9, 30, 66],
[8, 31, 58],
[8, 32, 57],
[8, 33, 56],
@ -1274,7 +1274,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[8, 40, 50],
[9, 41, 49],
[9, 42, 54],
[8, 43, 46],
[8, 43, 47],
[8, 44, 45],
[8, 48, 48],
[8, 50, 40],
@ -1287,9 +1287,9 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
[8, 71, 19],
[9, 72, 18],
[8, 79, 17],
[8, 81, 8],
[8, 81, 9],
[8, 89, 0],
],
], # 3
[
[7, 0, 28],
[8, 1, 27],
@ -1794,22 +1794,23 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
ir_weights_with_zeros = np.zeros(45, dtype=int)
ir_weights_with_zeros[: len(ir_weights)] = ir_weights
np.testing.assert_equal(ref_triplets[i][params[3]], ir_triplets)
np.testing.assert_equal(ref_ir_weights[i][params[3]], ir_weights)
if i == 1 and params[0]:
# print("{")
# for j, tp in enumerate(triplets_with_zeros):
# print("%d, %d, %d, " % tuple(tp), end="")
# if (j + 1) % 5 == 0:
# print("&")
print(", ".join(["%d" % x for x in ir_weights_with_zeros]))
# print("},")
# print(len(ir_triplets))
# print("[")
# print("[")
# for tp in ir_triplets:
# print("[%d, %d, %d]," % tuple(tp))
# print("],")
np.testing.assert_equal(ref_triplets[i][params[3]], ir_triplets)
np.testing.assert_equal(ref_ir_weights[i][params[3]], ir_weights)
# if i == 1 and params[0]:
# print("{")
# for j, tp in enumerate(triplets_with_zeros):
# print("%d, %d, %d, " % tuple(tp), end="")
# if (j + 1) % 5 == 0:
# print("&")
# print(", ".join(["%d" % x for x in ir_weights_with_zeros]))
# print("},")
# print(len(ir_triplets))
# print("[")
# print("]")
# print("[")
# print(",".join(["%d" % x for x in ir_weights]))