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-04-30 21:36:49 +09:00
parent 7806f7eb67 2e4c0b5d69
commit f0c26cec8a
52 changed files with 4551 additions and 4192 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.9.9
rev: v0.11.7
hooks:
- id: ruff
args: [ "--fix", "--show-fixes" ]

219
c/_phono3py.cpp
View File

@ -65,15 +65,15 @@ void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
Darray *fc3_normal_squared;
Darray *freqs;
_lapack_complex_double *eigvecs;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t num_triplets;
char *g_zero;
int64_t(*bz_grid_addresses)[3];
int64_t (*bz_grid_addresses)[3];
int64_t *D_diag;
int64_t(*Q)[3];
int64_t (*Q)[3];
double *fc3;
double(*svecs)[3];
int64_t(*multi)[2];
double (*svecs)[3];
int64_t (*multi)[2];
double *masses;
char *all_shortest;
int64_t *p2s;
@ -88,23 +88,23 @@ void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
/* npy_cdouble and lapack_complex_double may not be compatible. */
/* So eigenvectors should not be used in Python side */
eigvecs = (_lapack_complex_double *)py_eigenvectors.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
num_triplets = (int64_t)py_triplets.shape(0);
g_zero = (char *)py_g_zero.data();
bz_grid_addresses = (int64_t(*)[3])py_bz_grid_addresses.data();
bz_grid_addresses = (int64_t (*)[3])py_bz_grid_addresses.data();
D_diag = (int64_t *)py_D_diag.data();
Q = (int64_t(*)[3])py_Q.data();
Q = (int64_t (*)[3])py_Q.data();
fc3 = (double *)py_fc3.data();
if (py_fc3.shape(0) == py_fc3.shape(1)) {
is_compact_fc3 = 0;
} else {
is_compact_fc3 = 1;
}
svecs = (double(*)[3])py_svecs.data();
svecs = (double (*)[3])py_svecs.data();
for (i = 0; i < 2; i++) {
multi_dims[i] = py_multi.shape(i);
}
multi = (int64_t(*)[2])py_multi.data();
multi = (int64_t (*)[2])py_multi.data();
masses = (double *)py_masses.data();
p2s = (int64_t *)py_p2s_map.data();
s2p = (int64_t *)py_s2p_map.data();
@ -131,74 +131,74 @@ void py_get_pp_collision(
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_band_indices, nb::ndarray<> py_temperatures, int64_t is_NU,
int64_t symmetrize_fc3_q, int64_t make_r0_average,
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,
int64_t openmp_per_triplets) {
double *gamma;
int64_t(*relative_grid_address)[4][3];
int64_t (*relative_grid_address)[4][3];
double *frequencies;
_lapack_complex_double *eigenvectors;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t num_triplets;
int64_t *triplet_weights;
int64_t(*bz_grid_addresses)[3];
int64_t (*bz_grid_addresses)[3];
int64_t *bz_map;
int64_t *D_diag;
int64_t(*Q)[3];
int64_t (*Q)[3];
double *fc3;
double(*svecs)[3];
int64_t(*multi)[2];
double (*svecs)[3];
int64_t (*multi)[2];
double *masses;
int64_t *p2s;
int64_t *s2p;
Larray *band_indices;
Darray *temperatures;
Darray *temperatures_THz;
char *all_shortest;
int64_t multi_dims[2];
int64_t i;
int64_t is_compact_fc3;
gamma = (double *)py_gamma.data();
relative_grid_address = (int64_t(*)[4][3])py_relative_grid_address.data();
relative_grid_address = (int64_t (*)[4][3])py_relative_grid_address.data();
frequencies = (double *)py_frequencies.data();
eigenvectors = (_lapack_complex_double *)py_eigenvectors.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
num_triplets = (int64_t)py_triplets.shape(0);
triplet_weights = (int64_t *)py_triplet_weights.data();
bz_grid_addresses = (int64_t(*)[3])py_bz_grid_addresses.data();
bz_grid_addresses = (int64_t (*)[3])py_bz_grid_addresses.data();
bz_map = (int64_t *)py_bz_map.data();
D_diag = (int64_t *)py_D_diag.data();
Q = (int64_t(*)[3])py_Q.data();
Q = (int64_t (*)[3])py_Q.data();
fc3 = (double *)py_fc3.data();
if (py_fc3.shape(0) == py_fc3.shape(1)) {
is_compact_fc3 = 0;
} else {
is_compact_fc3 = 1;
}
svecs = (double(*)[3])py_svecs.data();
svecs = (double (*)[3])py_svecs.data();
for (i = 0; i < 2; i++) {
multi_dims[i] = py_multi.shape(i);
}
multi = (int64_t(*)[2])py_multi.data();
multi = (int64_t (*)[2])py_multi.data();
masses = (double *)py_masses.data();
p2s = (int64_t *)py_p2s_map.data();
s2p = (int64_t *)py_s2p_map.data();
band_indices = convert_to_larray(py_band_indices);
temperatures = convert_to_darray(py_temperatures);
temperatures_THz = convert_to_darray(py_temperatures_THz);
all_shortest = (char *)py_all_shortest.data();
ph3py_get_pp_collision(
gamma, relative_grid_address, frequencies, eigenvectors, triplets,
num_triplets, triplet_weights, bz_grid_addresses, bz_map, bz_grid_type,
D_diag, Q, fc3, is_compact_fc3, svecs, multi_dims, multi, masses, p2s,
s2p, band_indices, temperatures, is_NU, symmetrize_fc3_q,
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;
free(temperatures);
temperatures = NULL;
free(temperatures_THz);
temperatures_THz = NULL;
}
void py_get_pp_collision_with_sigma(
@ -209,26 +209,26 @@ void py_get_pp_collision_with_sigma(
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,
nb::ndarray<> py_temperatures, int64_t is_NU, int64_t symmetrize_fc3_q,
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, int64_t openmp_per_triplets) {
double *gamma;
double *frequencies;
_lapack_complex_double *eigenvectors;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t num_triplets;
int64_t *triplet_weights;
int64_t(*bz_grid_addresses)[3];
int64_t (*bz_grid_addresses)[3];
int64_t *D_diag;
int64_t(*Q)[3];
int64_t (*Q)[3];
double *fc3;
double(*svecs)[3];
int64_t(*multi)[2];
double (*svecs)[3];
int64_t (*multi)[2];
double *masses;
int64_t *p2s;
int64_t *s2p;
Larray *band_indices;
Darray *temperatures;
Darray *temperatures_THz;
char *all_shortest;
int64_t multi_dims[2];
int64_t i;
@ -237,47 +237,47 @@ void py_get_pp_collision_with_sigma(
gamma = (double *)py_gamma.data();
frequencies = (double *)py_frequencies.data();
eigenvectors = (_lapack_complex_double *)py_eigenvectors.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
num_triplets = (int64_t)py_triplets.shape(0);
triplet_weights = (int64_t *)py_triplet_weights.data();
bz_grid_addresses = (int64_t(*)[3])py_bz_grid_addresses.data();
bz_grid_addresses = (int64_t (*)[3])py_bz_grid_addresses.data();
D_diag = (int64_t *)py_D_diag.data();
Q = (int64_t(*)[3])py_Q.data();
Q = (int64_t (*)[3])py_Q.data();
fc3 = (double *)py_fc3.data();
if (py_fc3.shape(0) == py_fc3.shape(1)) {
is_compact_fc3 = 0;
} else {
is_compact_fc3 = 1;
}
svecs = (double(*)[3])py_svecs.data();
svecs = (double (*)[3])py_svecs.data();
for (i = 0; i < 2; i++) {
multi_dims[i] = py_multi.shape(i);
}
multi = (int64_t(*)[2])py_multi.data();
multi = (int64_t (*)[2])py_multi.data();
masses = (double *)py_masses.data();
p2s = (int64_t *)py_p2s_map.data();
s2p = (int64_t *)py_s2p_map.data();
band_indices = convert_to_larray(py_band_indices);
temperatures = convert_to_darray(py_temperatures);
temperatures_THz = convert_to_darray(py_temperatures_THz);
all_shortest = (char *)py_all_shortest.data();
ph3py_get_pp_collision_with_sigma(
gamma, sigma, sigma_cutoff, frequencies, eigenvectors, triplets,
num_triplets, triplet_weights, bz_grid_addresses, D_diag, Q, fc3,
is_compact_fc3, svecs, multi_dims, multi, masses, p2s, s2p,
band_indices, temperatures, is_NU, symmetrize_fc3_q, make_r0_average,
all_shortest, cutoff_frequency, openmp_per_triplets);
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;
free(temperatures);
temperatures = NULL;
free(temperatures_THz);
temperatures_THz = NULL;
}
void py_get_imag_self_energy_with_g(
nb::ndarray<> py_gamma, nb::ndarray<> py_fc3_normal_squared,
nb::ndarray<> py_triplets, nb::ndarray<> py_triplet_weights,
nb::ndarray<> py_frequencies, double temperature, nb::ndarray<> py_g,
nb::ndarray<> py_frequencies, double temperature_THz, nb::ndarray<> py_g,
nb::ndarray<> py_g_zero, double cutoff_frequency,
int64_t frequency_point_index) {
Darray *fc3_normal_squared;
@ -285,7 +285,7 @@ void py_get_imag_self_energy_with_g(
double *g;
char *g_zero;
double *frequencies;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t *triplet_weights;
int64_t num_frequency_points;
@ -294,13 +294,13 @@ void py_get_imag_self_energy_with_g(
g = (double *)py_g.data();
g_zero = (char *)py_g_zero.data();
frequencies = (double *)py_frequencies.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
triplet_weights = (int64_t *)py_triplet_weights.data();
num_frequency_points = (int64_t)py_g.shape(2);
ph3py_get_imag_self_energy_at_bands_with_g(
gamma, fc3_normal_squared, frequencies, triplets, triplet_weights, g,
g_zero, temperature, cutoff_frequency, num_frequency_points,
g_zero, temperature_THz, cutoff_frequency, num_frequency_points,
frequency_point_index);
free(fc3_normal_squared);
@ -312,7 +312,7 @@ void py_get_detailed_imag_self_energy_with_g(
nb::ndarray<> py_gamma_U, nb::ndarray<> py_fc3_normal_squared,
nb::ndarray<> py_triplets, nb::ndarray<> py_triplet_weights,
nb::ndarray<> py_bz_grid_addresses, nb::ndarray<> py_frequencies,
double temperature, nb::ndarray<> py_g, nb::ndarray<> py_g_zero,
double temperature_THz, nb::ndarray<> py_g, nb::ndarray<> py_g_zero,
double cutoff_frequency) {
Darray *fc3_normal_squared;
double *gamma_detail;
@ -321,9 +321,9 @@ void py_get_detailed_imag_self_energy_with_g(
double *g;
char *g_zero;
double *frequencies;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t *triplet_weights;
int64_t(*bz_grid_addresses)[3];
int64_t (*bz_grid_addresses)[3];
fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
gamma_detail = (double *)py_gamma_detail.data();
@ -332,14 +332,14 @@ void py_get_detailed_imag_self_energy_with_g(
g = (double *)py_g.data();
g_zero = (char *)py_g_zero.data();
frequencies = (double *)py_frequencies.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
triplet_weights = (int64_t *)py_triplet_weights.data();
bz_grid_addresses = (int64_t(*)[3])py_bz_grid_addresses.data();
bz_grid_addresses = (int64_t (*)[3])py_bz_grid_addresses.data();
ph3py_get_detailed_imag_self_energy_at_bands_with_g(
gamma_detail, gamma_N, gamma_U, fc3_normal_squared, frequencies,
triplets, triplet_weights, bz_grid_addresses, g, g_zero, temperature,
cutoff_frequency);
triplets, triplet_weights, bz_grid_addresses, g, g_zero,
temperature_THz, cutoff_frequency);
free(fc3_normal_squared);
fc3_normal_squared = NULL;
@ -349,25 +349,25 @@ void py_get_real_self_energy_at_bands(
nb::ndarray<> py_shift, nb::ndarray<> py_fc3_normal_squared,
nb::ndarray<> py_triplets, nb::ndarray<> py_triplet_weights,
nb::ndarray<> py_frequencies, nb::ndarray<> py_band_indices,
double temperature, double epsilon, double unit_conversion_factor,
double temperature_THz, double epsilon, double unit_conversion_factor,
double cutoff_frequency) {
Darray *fc3_normal_squared;
double *shift;
double *frequencies;
int64_t *band_indices;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t *triplet_weights;
fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
shift = (double *)py_shift.data();
frequencies = (double *)py_frequencies.data();
band_indices = (int64_t *)py_band_indices.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
triplet_weights = (int64_t *)py_triplet_weights.data();
ph3py_get_real_self_energy_at_bands(
shift, fc3_normal_squared, band_indices, frequencies, triplets,
triplet_weights, epsilon, temperature, unit_conversion_factor,
triplet_weights, epsilon, temperature_THz, unit_conversion_factor,
cutoff_frequency);
free(fc3_normal_squared);
@ -378,26 +378,26 @@ void py_get_real_self_energy_at_frequency_point(
nb::ndarray<> py_shift, double frequency_point,
nb::ndarray<> py_fc3_normal_squared, nb::ndarray<> py_triplets,
nb::ndarray<> py_triplet_weights, nb::ndarray<> py_frequencies,
nb::ndarray<> py_band_indices, double temperature, double epsilon,
nb::ndarray<> py_band_indices, double temperature_THz, double epsilon,
double unit_conversion_factor, double cutoff_frequency) {
Darray *fc3_normal_squared;
double *shift;
double *frequencies;
int64_t *band_indices;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t *triplet_weights;
fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
shift = (double *)py_shift.data();
frequencies = (double *)py_frequencies.data();
band_indices = (int64_t *)py_band_indices.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
triplet_weights = (int64_t *)py_triplet_weights.data();
ph3py_get_real_self_energy_at_frequency_point(
shift, frequency_point, fc3_normal_squared, band_indices, frequencies,
triplets, triplet_weights, epsilon, temperature, unit_conversion_factor,
cutoff_frequency);
triplets, triplet_weights, epsilon, temperature_THz,
unit_conversion_factor, cutoff_frequency);
free(fc3_normal_squared);
fc3_normal_squared = NULL;
@ -408,13 +408,13 @@ void py_get_collision_matrix(
nb::ndarray<> py_frequencies, nb::ndarray<> py_g, nb::ndarray<> py_triplets,
nb::ndarray<> py_triplets_map, nb::ndarray<> py_map_q,
nb::ndarray<> py_rotated_grid_points, nb::ndarray<> py_rotations_cartesian,
double temperature, double unit_conversion_factor,
double temperature_THz, double unit_conversion_factor,
double cutoff_frequency) {
Darray *fc3_normal_squared;
double *collision_matrix;
double *g;
double *frequencies;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t *triplets_map;
int64_t *map_q;
int64_t *rotated_grid_points;
@ -425,7 +425,7 @@ void py_get_collision_matrix(
collision_matrix = (double *)py_collision_matrix.data();
g = (double *)py_g.data();
frequencies = (double *)py_frequencies.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
triplets_map = (int64_t *)py_triplets_map.data();
num_gp = (int64_t)py_triplets_map.shape(0);
map_q = (int64_t *)py_map_q.data();
@ -440,7 +440,7 @@ void py_get_collision_matrix(
ph3py_get_collision_matrix(collision_matrix, fc3_normal_squared,
frequencies, triplets, triplets_map, map_q,
rotated_grid_points, rotations_cartesian, g,
num_ir_gp, num_gp, num_rot, temperature,
num_ir_gp, num_gp, num_rot, temperature_THz,
unit_conversion_factor, cutoff_frequency);
free(fc3_normal_squared);
@ -450,13 +450,14 @@ void py_get_collision_matrix(
void py_get_reducible_collision_matrix(
nb::ndarray<> py_collision_matrix, nb::ndarray<> py_fc3_normal_squared,
nb::ndarray<> py_frequencies, nb::ndarray<> py_g, nb::ndarray<> py_triplets,
nb::ndarray<> py_triplets_map, nb::ndarray<> py_map_q, double temperature,
double unit_conversion_factor, double cutoff_frequency) {
nb::ndarray<> py_triplets_map, nb::ndarray<> py_map_q,
double temperature_THz, double unit_conversion_factor,
double cutoff_frequency) {
Darray *fc3_normal_squared;
double *collision_matrix;
double *g;
double *frequencies;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
int64_t *triplets_map;
int64_t num_gp;
int64_t *map_q;
@ -465,14 +466,14 @@ void py_get_reducible_collision_matrix(
collision_matrix = (double *)py_collision_matrix.data();
g = (double *)py_g.data();
frequencies = (double *)py_frequencies.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
triplets_map = (int64_t *)py_triplets_map.data();
num_gp = (int64_t)py_triplets_map.shape(0);
map_q = (int64_t *)py_map_q.data();
ph3py_get_reducible_collision_matrix(
collision_matrix, fc3_normal_squared, frequencies, triplets,
triplets_map, map_q, g, num_gp, temperature, unit_conversion_factor,
triplets_map, map_q, g, num_gp, temperature_THz, unit_conversion_factor,
cutoff_frequency);
free(fc3_normal_squared);
@ -544,21 +545,21 @@ void py_rotate_delta_fc2s(nb::ndarray<> py_fc3, nb::ndarray<> py_delta_fc2s,
nb::ndarray<> py_inv_U,
nb::ndarray<> py_site_sym_cart,
nb::ndarray<> py_rot_map_syms) {
double(*fc3)[3][3][3];
double(*delta_fc2s)[3][3];
double (*fc3)[3][3][3];
double (*delta_fc2s)[3][3];
double *inv_U;
double(*site_sym_cart)[3][3];
double (*site_sym_cart)[3][3];
int64_t *rot_map_syms;
int64_t num_atom, num_disp, num_site_sym;
/* (num_atom, num_atom, 3, 3, 3) */
fc3 = (double(*)[3][3][3])py_fc3.data();
fc3 = (double (*)[3][3][3])py_fc3.data();
/* (n_u1, num_atom, num_atom, 3, 3) */
delta_fc2s = (double(*)[3][3])py_delta_fc2s.data();
delta_fc2s = (double (*)[3][3])py_delta_fc2s.data();
/* (3, n_u1 * n_sym) */
inv_U = (double *)py_inv_U.data();
/* (n_sym, 3, 3) */
site_sym_cart = (double(*)[3][3])py_site_sym_cart.data();
site_sym_cart = (double (*)[3][3])py_site_sym_cart.data();
/* (n_sym, natom) */
rot_map_syms = (int64_t *)py_rot_map_syms.data();
@ -695,11 +696,11 @@ void py_transpose_compact_fc3(nb::ndarray<> py_fc3,
void py_get_thm_relative_grid_address(nb::ndarray<> py_relative_grid_address,
nb::ndarray<> py_reciprocal_lattice_py) {
int64_t(*relative_grid_address)[4][3];
double(*reciprocal_lattice)[3];
int64_t (*relative_grid_address)[4][3];
double (*reciprocal_lattice)[3];
relative_grid_address = (int64_t(*)[4][3])py_relative_grid_address.data();
reciprocal_lattice = (double(*)[3])py_reciprocal_lattice_py.data();
relative_grid_address = (int64_t (*)[4][3])py_relative_grid_address.data();
reciprocal_lattice = (double (*)[3])py_reciprocal_lattice_py.data();
ph3py_get_relative_grid_address(relative_grid_address, reciprocal_lattice);
}
@ -714,18 +715,18 @@ void py_get_neighboring_grid_points(nb::ndarray<> py_relative_grid_points,
int64_t *relative_grid_points;
int64_t *grid_points;
int64_t num_grid_points, num_relative_grid_address;
int64_t(*relative_grid_address)[3];
int64_t (*relative_grid_address)[3];
int64_t *D_diag;
int64_t(*bz_grid_address)[3];
int64_t (*bz_grid_address)[3];
int64_t *bz_map;
relative_grid_points = (int64_t *)py_relative_grid_points.data();
grid_points = (int64_t *)py_grid_points.data();
num_grid_points = (int64_t)py_grid_points.shape(0);
relative_grid_address = (int64_t(*)[3])py_relative_grid_address.data();
relative_grid_address = (int64_t (*)[3])py_relative_grid_address.data();
num_relative_grid_address = (int64_t)py_relative_grid_address.shape(0);
D_diag = (int64_t *)py_D_diag.data();
bz_grid_address = (int64_t(*)[3])py_bz_grid_address.data();
bz_grid_address = (int64_t (*)[3])py_bz_grid_address.data();
bz_map = (int64_t *)py_bz_map.data();
ph3py_get_neighboring_gird_points(
@ -743,10 +744,10 @@ void py_get_thm_integration_weights_at_grid_points(
double *iw;
double *frequency_points;
int64_t num_frequency_points, num_band, num_gp;
int64_t(*relative_grid_address)[4][3];
int64_t (*relative_grid_address)[4][3];
int64_t *D_diag;
int64_t *grid_points;
int64_t(*bz_grid_address)[3];
int64_t (*bz_grid_address)[3];
int64_t *bz_map;
int64_t *gp2irgp_map;
double *frequencies;
@ -754,11 +755,11 @@ void py_get_thm_integration_weights_at_grid_points(
iw = (double *)py_iw.data();
frequency_points = (double *)py_frequency_points.data();
num_frequency_points = (int64_t)py_frequency_points.shape(0);
relative_grid_address = (int64_t(*)[4][3])py_relative_grid_address.data();
relative_grid_address = (int64_t (*)[4][3])py_relative_grid_address.data();
D_diag = (int64_t *)py_D_diag.data();
grid_points = (int64_t *)py_grid_points.data();
num_gp = (int64_t)py_grid_points.shape(0);
bz_grid_address = (int64_t(*)[3])py_bz_grid_address.data();
bz_grid_address = (int64_t (*)[3])py_bz_grid_address.data();
bz_map = (int64_t *)py_bz_map.data();
gp2irgp_map = (int64_t *)py_gp2irgp_map.data();
frequencies = (double *)py_frequencies.data();
@ -777,14 +778,14 @@ int64_t py_tpl_get_triplets_reciprocal_mesh_at_q(
int64_t *map_triplets;
int64_t *map_q;
int64_t *D_diag;
int64_t(*rot)[3][3];
int64_t (*rot)[3][3];
int64_t num_rot;
int64_t num_ir;
map_triplets = (int64_t *)py_map_triplets.data();
map_q = (int64_t *)py_map_q.data();
D_diag = (int64_t *)py_D_diag.data();
rot = (int64_t(*)[3][3])py_rotations.data();
rot = (int64_t (*)[3][3])py_rotations.data();
num_rot = (int64_t)py_rotations.shape(0);
num_ir = ph3py_get_triplets_reciprocal_mesh_at_q(
@ -801,22 +802,22 @@ int64_t py_tpl_get_BZ_triplets_at_q(nb::ndarray<> py_triplets,
nb::ndarray<> py_map_triplets,
nb::ndarray<> py_D_diag, nb::ndarray<> py_Q,
int64_t bz_grid_type) {
int64_t(*triplets)[3];
int64_t(*bz_grid_address)[3];
int64_t (*triplets)[3];
int64_t (*bz_grid_address)[3];
int64_t *bz_map;
int64_t *map_triplets;
int64_t num_map_triplets;
int64_t *D_diag;
int64_t(*Q)[3];
int64_t (*Q)[3];
int64_t num_ir;
triplets = (int64_t(*)[3])py_triplets.data();
bz_grid_address = (int64_t(*)[3])py_bz_grid_address.data();
triplets = (int64_t (*)[3])py_triplets.data();
bz_grid_address = (int64_t (*)[3])py_bz_grid_address.data();
bz_map = (int64_t *)py_bz_map.data();
map_triplets = (int64_t *)py_map_triplets.data();
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();
Q = (int64_t (*)[3])py_Q.data();
num_ir = ph3py_get_BZ_triplets_at_q(triplets, grid_point, bz_grid_address,
bz_map, map_triplets, num_map_triplets,
@ -835,10 +836,10 @@ void py_get_triplets_integration_weights(
double *iw;
char *iw_zero;
double *frequency_points;
int64_t(*relative_grid_address)[4][3];
int64_t (*relative_grid_address)[4][3];
int64_t *D_diag;
int64_t(*triplets)[3];
int64_t(*bz_grid_addresses)[3];
int64_t (*triplets)[3];
int64_t (*bz_grid_addresses)[3];
int64_t *bz_map;
double *frequencies1, *frequencies2;
int64_t num_band0, num_band1, num_band2, num_triplets;
@ -847,11 +848,11 @@ void py_get_triplets_integration_weights(
iw_zero = (char *)py_iw_zero.data();
frequency_points = (double *)py_frequency_points.data();
num_band0 = (int64_t)py_frequency_points.shape(0);
relative_grid_address = (int64_t(*)[4][3])py_relative_grid_address.data();
relative_grid_address = (int64_t (*)[4][3])py_relative_grid_address.data();
D_diag = (int64_t *)py_D_diag.data();
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
num_triplets = (int64_t)py_triplets.shape(0);
bz_grid_addresses = (int64_t(*)[3])py_bz_grid_addresses.data();
bz_grid_addresses = (int64_t (*)[3])py_bz_grid_addresses.data();
bz_map = (int64_t *)py_bz_map.data();
frequencies1 = (double *)py_frequencies1.data();
frequencies2 = (double *)py_frequencies2.data();
@ -871,7 +872,7 @@ void py_get_triplets_integration_weights_with_sigma(
double *iw;
char *iw_zero;
double *frequency_points;
int64_t(*triplets)[3];
int64_t (*triplets)[3];
double *frequencies;
int64_t num_band0, num_band, num_iw, num_triplets;
@ -879,7 +880,7 @@ void py_get_triplets_integration_weights_with_sigma(
iw_zero = (char *)py_iw_zero.data();
frequency_points = (double *)py_frequency_points.data();
num_band0 = (int64_t)py_frequency_points.shape(0);
triplets = (int64_t(*)[3])py_triplets.data();
triplets = (int64_t (*)[3])py_triplets.data();
num_triplets = (int64_t)py_triplets.shape(0);
frequencies = (double *)py_frequencies.data();
num_band = (int64_t)py_frequencies.shape(1);

62
c/collision_matrix.c
View File

@ -49,7 +49,7 @@ static void get_collision_matrix(
const int64_t num_gp, const int64_t *map_q, const int64_t *rot_grid_points,
const int64_t num_ir_gp, const int64_t num_rot,
const double *rotations_cartesian, const double *g,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency);
static void get_collision_matrix_at_gp(
double *collision_matrix, const double *fc3_normal_squared,
@ -58,25 +58,25 @@ static void get_collision_matrix_at_gp(
const int64_t *map_q, const int64_t *rot_grid_points,
const int64_t num_ir_gp, const int64_t num_rot,
const double *rotations_cartesian, const double *g,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency, const int64_t *gp2tp_map, const int64_t i);
static void get_reducible_collision_matrix(
double *collision_matrix, const double *fc3_normal_squared,
const int64_t num_band0, const int64_t num_band, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplets_map,
const int64_t num_gp, const int64_t *map_q, const double *g,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency);
static void get_reducible_collision_matrix_at_gp(
double *collision_matrix, const double *fc3_normal_squared,
const int64_t num_band0, const int64_t num_band, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplets_map,
const int64_t num_gp, const int64_t *map_q, const double *g,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency, const int64_t *gp2tp_map, const int64_t i);
static int64_t get_inv_sinh(double *inv_sinh, const int64_t gp,
const double temperature, const double *frequencies,
const int64_t triplet[3],
const double temperature_THz,
const double *frequencies, const int64_t triplet[3],
const int64_t *triplets_map, const int64_t *map_q,
const int64_t num_band,
const double cutoff_frequency);
@ -89,7 +89,7 @@ void col_get_collision_matrix(
const int64_t *triplets_map, const int64_t *map_q,
const int64_t *rot_grid_points, const double *rotations_cartesian,
const double *g, const int64_t num_ir_gp, const int64_t num_gp,
const int64_t num_rot, const double temperature,
const int64_t num_rot, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency) {
int64_t num_triplets, num_band0, num_band;
@ -97,19 +97,19 @@ void col_get_collision_matrix(
num_band0 = fc3_normal_squared->dims[1];
num_band = fc3_normal_squared->dims[2];
get_collision_matrix(collision_matrix, fc3_normal_squared->data, num_band0,
num_band, frequencies, triplets, triplets_map, num_gp,
map_q, rot_grid_points, num_ir_gp, num_rot,
rotations_cartesian,
g + 2 * num_triplets * num_band0 * num_band * num_band,
temperature, unit_conversion_factor, cutoff_frequency);
get_collision_matrix(
collision_matrix, fc3_normal_squared->data, num_band0, num_band,
frequencies, triplets, triplets_map, num_gp, map_q, rot_grid_points,
num_ir_gp, num_rot, rotations_cartesian,
g + 2 * num_triplets * num_band0 * num_band * num_band, temperature_THz,
unit_conversion_factor, cutoff_frequency);
}
void col_get_reducible_collision_matrix(
double *collision_matrix, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplets_map, const int64_t *map_q, const double *g,
const int64_t num_gp, const double temperature,
const int64_t num_gp, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency) {
int64_t num_triplets, num_band, num_band0;
@ -120,7 +120,7 @@ void col_get_reducible_collision_matrix(
get_reducible_collision_matrix(
collision_matrix, fc3_normal_squared->data, num_band0, num_band,
frequencies, triplets, triplets_map, num_gp, map_q,
g + 2 * num_triplets * num_band0 * num_band * num_band, temperature,
g + 2 * num_triplets * num_band0 * num_band * num_band, temperature_THz,
unit_conversion_factor, cutoff_frequency);
}
@ -131,7 +131,7 @@ static void get_collision_matrix(
const int64_t num_gp, const int64_t *map_q, const int64_t *rot_grid_points,
const int64_t num_ir_gp, const int64_t num_rot,
const double *rotations_cartesian, const double *g,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency) {
int64_t i;
int64_t *gp2tp_map;
@ -145,7 +145,7 @@ static void get_collision_matrix(
get_collision_matrix_at_gp(
collision_matrix, fc3_normal_squared, num_band0, num_band,
frequencies, triplets, triplets_map, map_q, rot_grid_points,
num_ir_gp, num_rot, rotations_cartesian, g, temperature,
num_ir_gp, num_rot, rotations_cartesian, g, temperature_THz,
unit_conversion_factor, cutoff_frequency, gp2tp_map, i);
}
@ -160,7 +160,7 @@ static void get_collision_matrix_at_gp(
const int64_t *map_q, const int64_t *rot_grid_points,
const int64_t num_ir_gp, const int64_t num_rot,
const double *rotations_cartesian, const double *g,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency, const int64_t *gp2tp_map, const int64_t i) {
int64_t j, k, l, m, n, ti, r_gp, swapped;
@ -171,9 +171,9 @@ static void get_collision_matrix_at_gp(
for (j = 0; j < num_rot; j++) {
r_gp = rot_grid_points[i * num_rot + j];
ti = gp2tp_map[triplets_map[r_gp]];
swapped =
get_inv_sinh(inv_sinh, r_gp, temperature, frequencies, triplets[ti],
triplets_map, map_q, num_band, cutoff_frequency);
swapped = get_inv_sinh(inv_sinh, r_gp, temperature_THz, frequencies,
triplets[ti], triplets_map, map_q, num_band,
cutoff_frequency);
for (k = 0; k < num_band0; k++) {
for (l = 0; l < num_band; l++) {
@ -219,7 +219,7 @@ static void get_reducible_collision_matrix(
const int64_t num_band0, const int64_t num_band, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplets_map,
const int64_t num_gp, const int64_t *map_q, const double *g,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency) {
int64_t i;
int64_t *gp2tp_map;
@ -232,8 +232,9 @@ static void get_reducible_collision_matrix(
for (i = 0; i < num_gp; i++) {
get_reducible_collision_matrix_at_gp(
collision_matrix, fc3_normal_squared, num_band0, num_band,
frequencies, triplets, triplets_map, num_gp, map_q, g, temperature,
unit_conversion_factor, cutoff_frequency, gp2tp_map, i);
frequencies, triplets, triplets_map, num_gp, map_q, g,
temperature_THz, unit_conversion_factor, cutoff_frequency,
gp2tp_map, i);
}
free(gp2tp_map);
@ -245,7 +246,7 @@ static void get_reducible_collision_matrix_at_gp(
const int64_t num_band0, const int64_t num_band, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplets_map,
const int64_t num_gp, const int64_t *map_q, const double *g,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency, const int64_t *gp2tp_map, const int64_t i) {
int64_t j, k, l, ti, swapped;
double collision;
@ -253,8 +254,9 @@ static void get_reducible_collision_matrix_at_gp(
inv_sinh = (double *)malloc(sizeof(double) * num_band);
ti = gp2tp_map[triplets_map[i]];
swapped = get_inv_sinh(inv_sinh, i, temperature, frequencies, triplets[ti],
triplets_map, map_q, num_band, cutoff_frequency);
swapped =
get_inv_sinh(inv_sinh, i, temperature_THz, frequencies, triplets[ti],
triplets_map, map_q, num_band, cutoff_frequency);
for (j = 0; j < num_band0; j++) {
for (k = 0; k < num_band; k++) {
@ -288,8 +290,8 @@ static void get_reducible_collision_matrix_at_gp(
}
static int64_t get_inv_sinh(double *inv_sinh, const int64_t gp,
const double temperature, const double *frequencies,
const int64_t triplet[3],
const double temperature_THz,
const double *frequencies, const int64_t triplet[3],
const int64_t *triplets_map, const int64_t *map_q,
const int64_t num_band,
const double cutoff_frequency) {
@ -310,7 +312,7 @@ static int64_t get_inv_sinh(double *inv_sinh, const int64_t gp,
for (i = 0; i < num_band; i++) {
f = frequencies[gp2 * num_band + i];
if (f > cutoff_frequency) {
inv_sinh[i] = funcs_inv_sinh_occupation(f, temperature);
inv_sinh[i] = funcs_inv_sinh_occupation(f, temperature_THz);
} else {
inv_sinh[i] = 0;
}

4
c/collision_matrix.h
View File

@ -45,13 +45,13 @@ void col_get_collision_matrix(
const int64_t *triplets_map, const int64_t *map_q,
const int64_t *rot_grid_points, const double *rotations_cartesian,
const double *g, const int64_t num_ir_gp, const int64_t num_gp,
const int64_t num_rot, const double temperature,
const int64_t num_rot, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency);
void col_get_reducible_collision_matrix(
double *collision_matrix, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplets_map, const int64_t *map_q, const double *g,
const int64_t num_gp, const double temperature,
const int64_t num_gp, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency);
#endif

9
c/funcs.c
View File

@ -38,17 +38,16 @@
#include "phonoc_const.h"
#define THZTOEVPARKB 47.992398658977166
#define INVSQRT2PI 0.3989422804014327
double funcs_bose_einstein(const double x, const double t) {
return 1.0 / (exp(THZTOEVPARKB * x / t) - 1);
double funcs_bose_einstein(const double x, const double temperature_THz) {
return 1.0 / (exp(x / temperature_THz) - 1);
}
double funcs_gaussian(const double x, const double sigma) {
return INVSQRT2PI / sigma * exp(-x * x / 2 / sigma / sigma);
}
double funcs_inv_sinh_occupation(const double x, const double t) {
return 1.0 / sinh(x * THZTOEVPARKB / 2 / t);
double funcs_inv_sinh_occupation(const double x, const double temperature_THz) {
return 1.0 / sinh(x / 2 / temperature_THz);
}

36
c/imag_self_energy_with_g.c
View File

@ -53,7 +53,7 @@ static void detailed_imag_self_energy_at_triplet(
const int64_t num_band0, const int64_t num_band,
const double *fc3_normal_squared, const double *frequencies,
const int64_t triplet[3], const double *g1, const double *g2_3,
const char *g_zero, const double *temperatures, const int64_t num_temps,
const char *g_zero, const double *temperatures_THz, const int64_t num_temps,
const double cutoff_frequency);
static double collect_detailed_imag_self_energy(
double *imag_self_energy, const int64_t num_band,
@ -72,11 +72,11 @@ void ise_get_imag_self_energy_with_g(
double *imag_self_energy, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const double *g, const char *g_zero,
const double temperature, const double cutoff_frequency,
const double temperature_THz, const double cutoff_frequency,
const int64_t num_frequency_points, const int64_t frequency_point_index) {
int64_t i, j, num_triplets, num_band0, num_band, num_band_prod;
int64_t num_g_pos, g_index_dims, g_index_shift;
int64_t(*g_pos)[4];
int64_t (*g_pos)[4];
double *ise;
int64_t at_a_frequency_point;
@ -118,7 +118,7 @@ void ise_get_imag_self_energy_with_g(
* ise_set_g_pos works for frequency points as bands.
* set_g_pos_frequency_point works for frequency sampling mode.
*/
g_pos = (int64_t(*)[4])malloc(sizeof(int64_t[4]) * num_band_prod);
g_pos = (int64_t (*)[4])malloc(sizeof(int64_t[4]) * num_band_prod);
if (at_a_frequency_point) {
num_g_pos = set_g_pos_frequency_point(
g_pos, num_band0, num_band,
@ -134,7 +134,7 @@ void ise_get_imag_self_energy_with_g(
triplets[i], triplet_weights[i],
g + i * g_index_dims + g_index_shift,
g + (i + num_triplets) * g_index_dims + g_index_shift, g_pos,
num_g_pos, &temperature, 1, cutoff_frequency, 0,
num_g_pos, &temperature_THz, 1, cutoff_frequency, 0,
at_a_frequency_point);
free(g_pos);
@ -160,7 +160,7 @@ void ise_get_detailed_imag_self_energy_with_g(
double *imag_self_energy_U, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const int64_t (*bz_grid_addresses)[3],
const double *g, const char *g_zero, const double temperature,
const double *g, const char *g_zero, const double temperature_THz,
const double cutoff_frequency) {
double *ise;
int64_t i, j, num_triplets, num_band0, num_band, num_band_prod;
@ -187,7 +187,7 @@ void ise_get_detailed_imag_self_energy_with_g(
num_band0, num_band, fc3_normal_squared->data + i * num_band_prod,
frequencies, triplets[i], g + i * num_band_prod,
g + (i + num_triplets) * num_band_prod, g_zero + i * num_band_prod,
&temperature, 1, cutoff_frequency);
&temperature_THz, 1, cutoff_frequency);
}
is_N = (int64_t *)malloc(sizeof(int64_t) * num_triplets);
@ -224,7 +224,7 @@ void ise_imag_self_energy_at_triplet(
const double *fc3_normal_squared, const double *frequencies,
const int64_t triplet[3], const int64_t triplet_weight, const double *g1,
const double *g2_3, const int64_t (*g_pos)[4], const int64_t num_g_pos,
const double *temperatures, const int64_t num_temps,
const double *temperatures_THz, const int64_t num_temps,
const double cutoff_frequency, const int64_t openmp_per_triplets,
const int64_t at_a_frequency_point) {
int64_t i, j;
@ -237,7 +237,7 @@ void ise_imag_self_energy_at_triplet(
for (i = 0; i < num_temps; i++) {
set_occupations(n1 + i * num_band, n2 + i * num_band, num_band,
temperatures[i], triplet, frequencies,
temperatures_THz[i], triplet, frequencies,
cutoff_frequency);
}
@ -259,7 +259,7 @@ void ise_imag_self_energy_at_triplet(
n2[j * num_band + g_pos[i][2]] < 0) {
ise_at_g_pos[i * num_temps + j] = 0;
} else {
if (temperatures[j] > 0) {
if (temperatures_THz[j] > 0) {
ise_at_g_pos[i * num_temps + j] =
((n1[j * num_band + g_pos[i][1]] +
n2[j * num_band + g_pos[i][2]] + 1) *
@ -359,7 +359,7 @@ static void detailed_imag_self_energy_at_triplet(
const int64_t num_band0, const int64_t num_band,
const double *fc3_normal_squared, const double *frequencies,
const int64_t triplet[3], const double *g1, const double *g2_3,
const char *g_zero, const double *temperatures, const int64_t num_temps,
const char *g_zero, const double *temperatures_THz, const int64_t num_temps,
const double cutoff_frequency) {
int64_t i, j, adrs_shift;
double *n1, *n2;
@ -371,12 +371,12 @@ static void detailed_imag_self_energy_at_triplet(
n2 = (double *)malloc(sizeof(double) * num_band);
for (i = 0; i < num_temps; i++) {
set_occupations(n1, n2, num_band, temperatures[i], triplet, frequencies,
cutoff_frequency);
set_occupations(n1, n2, num_band, temperatures_THz[i], triplet,
frequencies, cutoff_frequency);
for (j = 0; j < num_band0; j++) {
adrs_shift = j * num_band * num_band;
if (temperatures[i] > 0) {
if (temperatures_THz[i] > 0) {
imag_self_energy[i * num_band0 + j] =
collect_detailed_imag_self_energy(
detailed_imag_self_energy + adrs_shift, num_band,
@ -452,8 +452,8 @@ static double collect_detailed_imag_self_energy_0K(
}
static void set_occupations(double *n1, double *n2, const int64_t num_band,
const double temperature, const int64_t triplet[3],
const double *frequencies,
const double temperature_THz,
const int64_t triplet[3], const double *frequencies,
const double cutoff_frequency) {
int64_t j;
double f1, f2;
@ -462,12 +462,12 @@ static void set_occupations(double *n1, double *n2, const int64_t num_band,
f1 = frequencies[triplet[1] * num_band + j];
f2 = frequencies[triplet[2] * num_band + j];
if (f1 > cutoff_frequency) {
n1[j] = funcs_bose_einstein(f1, temperature);
n1[j] = funcs_bose_einstein(f1, temperature_THz);
} else {
n1[j] = -1;
}
if (f2 > cutoff_frequency) {
n2[j] = funcs_bose_einstein(f2, temperature);
n2[j] = funcs_bose_einstein(f2, temperature_THz);
} else {
n2[j] = -1;
}

6
c/imag_self_energy_with_g.h
View File

@ -44,21 +44,21 @@ void ise_get_imag_self_energy_with_g(
double *imag_self_energy, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const double *g, const char *g_zero,
const double temperature, const double cutoff_frequency,
const double temperature_THz, const double cutoff_frequency,
const int64_t num_frequency_points, const int64_t frequency_point_index);
void ise_get_detailed_imag_self_energy_with_g(
double *detailed_imag_self_energy, double *imag_self_energy_N,
double *imag_self_energy_U, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const int64_t (*bz_grid_addresses)[3],
const double *g, const char *g_zero, const double temperature,
const double *g, const char *g_zero, const double temperature_THz,
const double cutoff_frequency);
void ise_imag_self_energy_at_triplet(
double *imag_self_energy, const int64_t num_band0, const int64_t num_band,
const double *fc3_normal_squared, const double *frequencies,
const int64_t triplet[3], const int64_t triplet_weight, const double *g1,
const double *g2_3, const int64_t (*g_pos)[4], const int64_t num_g_pos,
const double *temperatures, const int64_t num_temps,
const double *temperatures_THz, const int64_t num_temps,
const double cutoff_frequency, const int64_t openmp_possible,
const int64_t at_a_frequency_point);
int64_t ise_set_g_pos(int64_t (*g_pos)[4], const int64_t num_band0,

46
c/phono3py.c
View File

@ -130,10 +130,10 @@ int64_t ph3py_get_pp_collision(
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, 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 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;
@ -174,7 +174,7 @@ int64_t ph3py_get_pp_collision(
(lapack_complex_double *)eigenvectors, triplets,
num_triplets, triplet_weights, bzgrid, fc3,
is_compact_fc3, atom_triplets, masses, band_indices,
temperatures, is_NU, symmetrize_fc3_q,
temperatures_THz, is_NU, symmetrize_fc3_q,
cutoff_frequency, openmp_per_triplets);
free(atom_triplets);
@ -195,10 +195,10 @@ int64_t ph3py_get_pp_collision_with_sigma(
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, 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 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;
@ -237,8 +237,8 @@ int64_t ph3py_get_pp_collision_with_sigma(
imag_self_energy, sigma, sigma_cutoff, frequencies,
(lapack_complex_double *)eigenvectors, triplets, num_triplets,
triplet_weights, bzgrid, fc3, is_compact_fc3, atom_triplets, masses,
band_indices, temperatures, is_NU, symmetrize_fc3_q, cutoff_frequency,
openmp_per_triplets);
band_indices, temperatures_THz, is_NU, symmetrize_fc3_q,
cutoff_frequency, openmp_per_triplets);
free(atom_triplets);
atom_triplets = NULL;
@ -253,11 +253,11 @@ 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],
const int64_t *triplet_weights, const double *g, const char *g_zero,
const double temperature, const double cutoff_frequency,
const double temperature_THz, const double cutoff_frequency,
const int64_t num_frequency_points, const int64_t frequency_point_index) {
ise_get_imag_self_energy_with_g(
imag_self_energy, fc3_normal_squared, frequencies, triplets,
triplet_weights, g, g_zero, temperature, cutoff_frequency,
triplet_weights, g, g_zero, temperature_THz, cutoff_frequency,
num_frequency_points, frequency_point_index);
}
@ -266,23 +266,23 @@ void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
double *imag_self_energy_U, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const int64_t (*bz_grid_addresses)[3],
const double *g, const char *g_zero, const double temperature,
const double *g, const char *g_zero, const double temperature_THz,
const double cutoff_frequency) {
ise_get_detailed_imag_self_energy_with_g(
detailed_imag_self_energy, imag_self_energy_N, imag_self_energy_U,
fc3_normal_squared, frequencies, triplets, triplet_weights,
bz_grid_addresses, g, g_zero, temperature, cutoff_frequency);
bz_grid_addresses, g, g_zero, temperature_THz, cutoff_frequency);
}
void ph3py_get_real_self_energy_at_bands(
double *real_self_energy, const Darray *fc3_normal_squared,
const int64_t *band_indices, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplet_weights,
const double epsilon, const double temperature,
const double epsilon, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency) {
rse_get_real_self_energy_at_bands(real_self_energy, fc3_normal_squared,
band_indices, frequencies, triplets,
triplet_weights, epsilon, temperature,
triplet_weights, epsilon, temperature_THz,
unit_conversion_factor, cutoff_frequency);
}
@ -291,11 +291,11 @@ void ph3py_get_real_self_energy_at_frequency_point(
const Darray *fc3_normal_squared, const int64_t *band_indices,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const double epsilon,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency) {
rse_get_real_self_energy_at_frequency_point(
real_self_energy, frequency_point, fc3_normal_squared, band_indices,
frequencies, triplets, triplet_weights, epsilon, temperature,
frequencies, triplets, triplet_weights, epsilon, temperature_THz,
unit_conversion_factor, cutoff_frequency);
}
@ -305,12 +305,12 @@ void ph3py_get_collision_matrix(
const int64_t *triplets_map, const int64_t *map_q,
const int64_t *rotated_grid_points, const double *rotations_cartesian,
const double *g, const int64_t num_ir_gp, const int64_t num_gp,
const int64_t num_rot, const double temperature,
const int64_t num_rot, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency) {
col_get_collision_matrix(collision_matrix, fc3_normal_squared, frequencies,
triplets, triplets_map, map_q, rotated_grid_points,
rotations_cartesian, g, num_ir_gp, num_gp, num_rot,
temperature, unit_conversion_factor,
temperature_THz, unit_conversion_factor,
cutoff_frequency);
}
@ -318,11 +318,11 @@ void ph3py_get_reducible_collision_matrix(
double *collision_matrix, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplets_map, const int64_t *map_q, const double *g,
const int64_t num_gp, const double temperature,
const int64_t num_gp, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency) {
col_get_reducible_collision_matrix(
collision_matrix, fc3_normal_squared, frequencies, triplets,
triplets_map, map_q, g, num_gp, temperature, unit_conversion_factor,
triplets_map, map_q, g, num_gp, temperature_THz, unit_conversion_factor,
cutoff_frequency);
}

28
c/phono3py.h
View File

@ -71,10 +71,10 @@ int64_t ph3py_get_pp_collision(
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, 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 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);
int64_t ph3py_get_pp_collision_with_sigma(
double *imag_self_energy, const double sigma, const double sigma_cutoff,
const double *frequencies, const _lapack_complex_double *eigenvectors,
@ -84,35 +84,35 @@ int64_t ph3py_get_pp_collision_with_sigma(
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, 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 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],
const int64_t *triplet_weights, const double *g, const char *g_zero,
const double temperature, const double cutoff_frequency,
const double temperature_THz, const double cutoff_frequency,
const int64_t num_frequency_points, const int64_t frequency_point_index);
void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
double *detailed_imag_self_energy, double *imag_self_energy_N,
double *imag_self_energy_U, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const int64_t (*bz_grid_addresses)[3],
const double *g, const char *g_zero, const double temperature,
const double *g, const char *g_zero, const double temperature_THz,
const double cutoff_frequency);
void ph3py_get_real_self_energy_at_bands(
double *real_self_energy, const Darray *fc3_normal_squared,
const int64_t *band_indices, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplet_weights,
const double epsilon, const double temperature,
const double epsilon, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency);
void ph3py_get_real_self_energy_at_frequency_point(
double *real_self_energy, const double frequency_point,
const Darray *fc3_normal_squared, const int64_t *band_indices,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const double epsilon,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency);
void ph3py_get_collision_matrix(
double *collision_matrix, const Darray *fc3_normal_squared,
@ -120,13 +120,13 @@ void ph3py_get_collision_matrix(
const int64_t *triplets_map, const int64_t *map_q,
const int64_t *rotated_grid_points, const double *rotations_cartesian,
const double *g, const int64_t num_ir_gp, const int64_t num_gp,
const int64_t num_rot, const double temperature,
const int64_t num_rot, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency);
void ph3py_get_reducible_collision_matrix(
double *collision_matrix, const Darray *fc3_normal_squared,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplets_map, const int64_t *map_q, const double *g,
const int64_t num_gp, const double temperature,
const int64_t num_gp, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency);
void ph3py_get_isotope_scattering_strength(
double *gamma, const int64_t grid_point, const int64_t *ir_grid_points,

1
c/phonoc_const.h
View File

@ -36,6 +36,5 @@
#define __phonoc_const_H__
#define M_2PI 6.283185307179586
#define KB 8.6173382568083159E-05
#endif

40
c/pp_collision.c
View File

@ -50,7 +50,7 @@
static void get_collision(
double *ise, const int64_t num_band0, const int64_t num_band,
const int64_t num_temps, const double *temperatures, const double *g,
const int64_t num_temps, const double *temperatures_THz, const double *g,
const char *g_zero, const double *frequencies,
const lapack_complex_double *eigenvectors, const int64_t triplet[3],
const int64_t triplet_weight, const RecgridConstBZGrid *bzgrid,
@ -72,9 +72,9 @@ void ppc_get_pp_collision(
const int64_t *triplet_weights, const RecgridConstBZGrid *bzgrid,
const double *fc3, const int64_t is_compact_fc3,
const AtomTriplets *atom_triplets, const double *masses,
const Larray *band_indices, const Darray *temperatures, const int64_t is_NU,
const int64_t symmetrize_fc3_q, const double cutoff_frequency,
const int64_t openmp_per_triplets) {
const Larray *band_indices, const Darray *temperatures_THz,
const int64_t is_NU, const int64_t symmetrize_fc3_q,
const double cutoff_frequency, const int64_t openmp_per_triplets) {
int64_t i;
int64_t num_band, num_band0, num_band_prod, num_temps;
double *ise, *freqs_at_gp, *g;
@ -89,7 +89,7 @@ void ppc_get_pp_collision(
num_band0 = band_indices->dims[0];
num_band = atom_triplets->multi_dims[1] * 3;
num_band_prod = num_band0 * num_band * num_band;
num_temps = temperatures->dims[0];
num_temps = temperatures_THz->dims[0];
ise =
(double *)malloc(sizeof(double) * num_triplets * num_temps * num_band0);
freqs_at_gp = (double *)malloc(sizeof(double) * num_band0);
@ -101,8 +101,8 @@ void ppc_get_pp_collision(
tpl_set_relative_grid_address(tp_relative_grid_address,
relative_grid_address, 2);
#ifdef _OPENMP
#pragma omp parallel for schedule(guided) private( \
g, g_zero) if (openmp_per_triplets)
#pragma omp parallel for schedule(guided) \
private(g, g_zero) if (openmp_per_triplets)
#endif
for (i = 0; i < num_triplets; i++) {
g = (double *)malloc(sizeof(double) * 2 * num_band_prod);
@ -115,7 +115,7 @@ void ppc_get_pp_collision(
num_band, 2, openmp_per_triplets);
get_collision(ise + i * num_temps * num_band0, num_band0, num_band,
num_temps, temperatures->data, g, g_zero, frequencies,
num_temps, temperatures_THz->data, g, g_zero, frequencies,
eigenvectors, triplets[i], triplet_weights[i], bzgrid,
fc3, is_compact_fc3, atom_triplets, masses,
band_indices->data, symmetrize_fc3_q, cutoff_frequency,
@ -143,9 +143,9 @@ void ppc_get_pp_collision_with_sigma(
const int64_t *triplet_weights, const RecgridConstBZGrid *bzgrid,
const double *fc3, const int64_t is_compact_fc3,
const AtomTriplets *atom_triplets, const double *masses,
const Larray *band_indices, const Darray *temperatures, const int64_t is_NU,
const int64_t symmetrize_fc3_q, const double cutoff_frequency,
const int64_t openmp_per_triplets) {
const Larray *band_indices, const Darray *temperatures_THz,
const int64_t is_NU, const int64_t symmetrize_fc3_q,
const double cutoff_frequency, const int64_t openmp_per_triplets) {
int64_t i;
int64_t num_band, num_band0, num_band_prod, num_temps;
int64_t const_adrs_shift;
@ -161,7 +161,7 @@ void ppc_get_pp_collision_with_sigma(
num_band0 = band_indices->dims[0];
num_band = atom_triplets->multi_dims[1] * 3;
num_band_prod = num_band0 * num_band * num_band;
num_temps = temperatures->dims[0];
num_temps = temperatures_THz->dims[0];
const_adrs_shift = num_band_prod;
ise =
@ -175,8 +175,8 @@ void ppc_get_pp_collision_with_sigma(
cutoff = sigma * sigma_cutoff;
#ifdef _OPENMP
#pragma omp parallel for schedule(guided) private( \
g, g_zero) if (openmp_per_triplets)
#pragma omp parallel for schedule(guided) \
private(g, g_zero) if (openmp_per_triplets)
#endif
for (i = 0; i < num_triplets; i++) {
g = (double *)malloc(sizeof(double) * 2 * num_band_prod);
@ -186,7 +186,7 @@ void ppc_get_pp_collision_with_sigma(
const_adrs_shift, frequencies, num_band, 2, 1);
get_collision(ise + i * num_temps * num_band0, num_band0, num_band,
num_temps, temperatures->data, g, g_zero, frequencies,
num_temps, temperatures_THz->data, g, g_zero, frequencies,
eigenvectors, triplets[i], triplet_weights[i], bzgrid,
fc3, is_compact_fc3, atom_triplets, masses,
band_indices->data, symmetrize_fc3_q, cutoff_frequency,
@ -209,7 +209,7 @@ void ppc_get_pp_collision_with_sigma(
static void get_collision(
double *ise, const int64_t num_band0, const int64_t num_band,
const int64_t num_temps, const double *temperatures, const double *g,
const int64_t num_temps, const double *temperatures_THz, const double *g,
const char *g_zero, const double *frequencies,
const lapack_complex_double *eigenvectors, const int64_t triplet[3],
const int64_t triplet_weight, const RecgridConstBZGrid *bzgrid,
@ -220,14 +220,14 @@ static void get_collision(
int64_t i;
int64_t num_band_prod, num_g_pos;
double *fc3_normal_squared;
int64_t(*g_pos)[4];
int64_t (*g_pos)[4];
fc3_normal_squared = NULL;
g_pos = NULL;
num_band_prod = num_band0 * num_band * num_band;
fc3_normal_squared = (double *)malloc(sizeof(double) * num_band_prod);
g_pos = (int64_t(*)[4])malloc(sizeof(int64_t[4]) * num_band_prod);
g_pos = (int64_t (*)[4])malloc(sizeof(int64_t[4]) * num_band_prod);
for (i = 0; i < num_band_prod; i++) {
fc3_normal_squared[i] = 0;
@ -243,8 +243,8 @@ static void get_collision(
ise_imag_self_energy_at_triplet(
ise, num_band0, num_band, fc3_normal_squared, frequencies, triplet,
triplet_weight, g, g + num_band_prod, g_pos, num_g_pos, temperatures,
num_temps, cutoff_frequency, openmp_per_triplets, 0);
triplet_weight, g, g + num_band_prod, g_pos, num_g_pos,
temperatures_THz, num_temps, cutoff_frequency, openmp_per_triplets, 0);
free(fc3_normal_squared);
fc3_normal_squared = NULL;

12
c/pp_collision.h
View File

@ -50,9 +50,9 @@ void ppc_get_pp_collision(
const int64_t *triplet_weights, const RecgridConstBZGrid *bzgrid,
const double *fc3, const int64_t is_compact_fc3,
const AtomTriplets *atom_triplets, const double *masses,
const Larray *band_indices, const Darray *temperatures, const int64_t is_NU,
const int64_t symmetrize_fc3_q, const double cutoff_frequency,
const int64_t openmp_per_triplets);
const Larray *band_indices, const Darray *temperatures_THz,
const int64_t is_NU, const int64_t symmetrize_fc3_q,
const double cutoff_frequency, const int64_t openmp_per_triplets);
void ppc_get_pp_collision_with_sigma(
double *imag_self_energy, const double sigma, const double sigma_cutoff,
const double *frequencies, const lapack_complex_double *eigenvectors,
@ -60,8 +60,8 @@ void ppc_get_pp_collision_with_sigma(
const int64_t *triplet_weights, const RecgridConstBZGrid *bzgrid,
const double *fc3, const int64_t is_compact_fc3,
const AtomTriplets *atom_triplets, const double *masses,
const Larray *band_indices, const Darray *temperatures, const int64_t is_NU,
const int64_t symmetrize_fc3_q, const double cutoff_frequency,
const int64_t openmp_per_triplets);
const Larray *band_indices, const Darray *temperatures_THz,
const int64_t is_NU, const int64_t symmetrize_fc3_q,
const double cutoff_frequency, const int64_t openmp_per_triplets);
#endif

28
c/real_self_energy.c
View File

@ -46,12 +46,12 @@ static double get_real_self_energy_at_band(
const int64_t band_index, const Darray *fc3_normal_squared,
const double fpoint, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplet_weights,
const double epsilon, const double temperature,
const double epsilon, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency);
static double sum_real_self_energy_at_band(
const int64_t num_band, const double *fc3_normal_squared,
const double fpoint, const double *freqs1, const double *freqs2,
const double epsilon, const double temperature,
const double epsilon, const double temperature_THz,
const double cutoff_frequency);
static double sum_real_self_energy_at_band_0K(
const int64_t num_band, const double *fc3_normal_squared,
@ -62,7 +62,7 @@ void rse_get_real_self_energy_at_bands(
double *real_self_energy, const Darray *fc3_normal_squared,
const int64_t *band_indices, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplet_weights,
const double epsilon, const double temperature,
const double epsilon, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency) {
int64_t i, num_band0, num_band, gp0;
double fpoint;
@ -79,8 +79,8 @@ void rse_get_real_self_energy_at_bands(
} else {
real_self_energy[i] = get_real_self_energy_at_band(
i, fc3_normal_squared, fpoint, frequencies, triplets,
triplet_weights, epsilon, temperature, unit_conversion_factor,
cutoff_frequency);
triplet_weights, epsilon, temperature_THz,
unit_conversion_factor, cutoff_frequency);
}
}
}
@ -90,7 +90,7 @@ void rse_get_real_self_energy_at_frequency_point(
const Darray *fc3_normal_squared, const int64_t *band_indices,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const double epsilon,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency) {
int64_t i, num_band0;
@ -103,8 +103,8 @@ void rse_get_real_self_energy_at_frequency_point(
} else {
real_self_energy[i] = get_real_self_energy_at_band(
i, fc3_normal_squared, frequency_point, frequencies, triplets,
triplet_weights, epsilon, temperature, unit_conversion_factor,
cutoff_frequency);
triplet_weights, epsilon, temperature_THz,
unit_conversion_factor, cutoff_frequency);
}
}
}
@ -113,7 +113,7 @@ static double get_real_self_energy_at_band(
const int64_t band_index, const Darray *fc3_normal_squared,
const double fpoint, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplet_weights,
const double epsilon, const double temperature,
const double epsilon, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency) {
int64_t i, num_triplets, num_band0, num_band, gp1, gp2;
double shift;
@ -129,14 +129,14 @@ static double get_real_self_energy_at_band(
for (i = 0; i < num_triplets; i++) {
gp1 = triplets[i][1];
gp2 = triplets[i][2];
if (temperature > 0) {
if (temperature_THz > 0) {
shift += sum_real_self_energy_at_band(
num_band,
fc3_normal_squared->data +
i * num_band0 * num_band * num_band +
band_index * num_band * num_band,
fpoint, frequencies + gp1 * num_band,
frequencies + gp2 * num_band, epsilon, temperature,
frequencies + gp2 * num_band, epsilon, temperature_THz,
cutoff_frequency) *
triplet_weights[i] * unit_conversion_factor;
} else {
@ -157,7 +157,7 @@ static double get_real_self_energy_at_band(
static double sum_real_self_energy_at_band(
const int64_t num_band, const double *fc3_normal_squared,
const double fpoint, const double *freqs1, const double *freqs2,
const double epsilon, const double temperature,
const double epsilon, const double temperature_THz,
const double cutoff_frequency) {
int64_t i, j;
double n1, n2, f1, f2, f3, f4, shift;
@ -166,10 +166,10 @@ static double sum_real_self_energy_at_band(
shift = 0;
for (i = 0; i < num_band; i++) {
if (freqs1[i] > cutoff_frequency) {
n1 = funcs_bose_einstein(freqs1[i], temperature);
n1 = funcs_bose_einstein(freqs1[i], temperature_THz);
for (j = 0; j < num_band; j++) {
if (freqs2[j] > cutoff_frequency) {
n2 = funcs_bose_einstein(freqs2[j], temperature);
n2 = funcs_bose_einstein(freqs2[j], temperature_THz);
f1 = fpoint + freqs1[i] + freqs2[j];
f2 = fpoint - freqs1[i] - freqs2[j];
f3 = fpoint - freqs1[i] + freqs2[j];

4
c/real_self_energy.h
View File

@ -45,12 +45,12 @@ void rse_get_real_self_energy_at_bands(
double *real_self_energy, const Darray *fc3_normal_squared,
const int64_t *band_indices, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplet_weights,
const double epsilon, const double temperature,
const double epsilon, const double temperature_THz,
const double unit_conversion_factor, const double cutoff_frequency);
void rse_get_real_self_energy_at_frequency_point(
double *real_self_energy, const double frequency_point,
const Darray *fc3_normal_squared, const int64_t *band_indices,
const double *frequencies, const int64_t (*triplets)[3],
const int64_t *triplet_weights, const double epsilon,
const double temperature, const double unit_conversion_factor,
const double temperature_THz, const double unit_conversion_factor,
const double cutoff_frequency);

4
doc/changelog.md
View File

@ -2,6 +2,10 @@
# Change Log
## Apr-30-2025: Version 3.15.1
- Release to follow the change of phonopy
## Mar-4-2025: Version 3.15.0
- Release to follow the change of phonopy

2
doc/conf.py
View File

@ -60,7 +60,7 @@ copyright = "2015, Atsushi Togo"
# The short X.Y version.
version = "3.15"
# The full version, including alpha/beta/rc tags.
release = "3.15.0"
release = "3.15.1"
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.

3
phono3py/api_isotope.py
View File

@ -35,7 +35,6 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.units import VaspToTHz
from phono3py.other.isotope import Isotope
@ -50,7 +49,7 @@ class Phono3pyIsotope:
mass_variances=None, # length of list is num_atom.
band_indices=None,
sigmas=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
use_grg=False,
symprec=1e-5,
cutoff_frequency=None,

9
phono3py/api_jointdos.py
View File

@ -35,9 +35,9 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import Primitive, Supercell
from phonopy.structure.symmetry import Symmetry
from phonopy.units import VaspToTHz
from phono3py.file_IO import write_joint_dos
from phono3py.phonon.grid import BZGrid
@ -65,7 +65,7 @@ class Phono3pyJointDos:
num_frequency_points=None,
num_points_in_batch=None,
temperatures=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
frequency_scale_factor=None,
use_grg=False,
SNF_coordinates="reciprocal",
@ -87,7 +87,10 @@ class Phono3pyJointDos:
else:
self._sigmas = sigmas
self._cutoff_frequency = cutoff_frequency
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().DefaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
self._is_mesh_symmetry = is_mesh_symmetry
self._is_symmetry = is_symmetry

13
phono3py/api_phono3py.py
View File

@ -58,6 +58,7 @@ from phonopy.interface.mlp import PhonopyMLP
from phonopy.interface.pypolymlp import (
PypolymlpParams,
)
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import (
Primitive,
@ -69,10 +70,9 @@ from phonopy.structure.cells import (
shape_supercell_matrix,
)
from phonopy.structure.symmetry import Symmetry
from phonopy.units import VaspToTHz
from phono3py.conductivity.direct_solution import get_thermal_conductivity_LBTE
from phono3py.conductivity.rta import get_thermal_conductivity_RTA
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.phono3py_yaml import Phono3pyYaml
from phono3py.phonon.grid import BZGrid
@ -149,7 +149,7 @@ class Phono3py:
primitive_matrix=None,
phonon_supercell_matrix=None,
cutoff_frequency=1e-4,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
is_symmetry=True,
is_mesh_symmetry=True,
use_grg=False,
@ -218,7 +218,10 @@ class Phono3py:
"""
self._symprec = symprec
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().DefaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._is_symmetry = is_symmetry
self._is_mesh_symmetry = is_mesh_symmetry
self._use_grg = use_grg

29
phono3py/conductivity/base.py
View File

@ -42,7 +42,7 @@ from typing import List, Optional, Tuple, Union
import numpy as np
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.phonon.thermal_properties import mode_cv
from phonopy.units import EV, Angstrom, Kb, THz, THzToEv
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
@ -50,9 +50,17 @@ from phono3py.phonon3.collision_matrix import CollisionMatrix
from phono3py.phonon3.imag_self_energy import ImagSelfEnergy
from phono3py.phonon3.interaction import Interaction
unit_to_WmK = (
(THz * Angstrom) ** 2 / (Angstrom**3) * EV / THz / (2 * np.pi)
) # 2pi comes from definition of lifetime.
def get_unit_to_WmK() -> float:
"""Return conversion factor to WmK."""
unit_to_WmK = (
(get_physical_units().THz * get_physical_units().Angstrom) ** 2
/ (get_physical_units().Angstrom ** 3)
* get_physical_units().EV
/ get_physical_units().THz
/ (2 * np.pi)
) # 2pi comes from definition of lifetime.
return unit_to_WmK
class HeatCapacityMixIn:
@ -96,15 +104,18 @@ class HeatCapacityMixIn:
"""
grid_point = self._grid_points[i_gp]
freqs = self._frequencies[grid_point][self._pp.band_indices] * THzToEv
cutoff = self._pp.cutoff_frequency * THzToEv
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 * Kb
condition = f < 100 * self._temperatures * get_physical_units().KB
cv[:, i] = np.where(
condition,
mode_cv(np.where(condition, self._temperatures, 10000), f),
@ -379,7 +390,7 @@ class ConductivityBase(ABC):
self._gamma_iso: Optional[np.ndarray] = None
volume = self._pp.primitive.volume
self._conversion_factor = unit_to_WmK / volume
self._conversion_factor = get_unit_to_WmK() / volume
self._averaged_pp_interaction = None
@ -903,7 +914,7 @@ class ConductivityBase(ABC):
for ll in range(num_band):
g_boundary[ll] = (
np.linalg.norm(gv[i_gp, ll])
* Angstrom
* get_physical_units().Angstrom
* 1e6
/ (4 * np.pi * self._boundary_mfp)
)

1743
phono3py/conductivity/direct_solution.py
View File

File diff suppressed because it is too large Load Diff

1352
phono3py/conductivity/direct_solution_base.py Normal file
View File

File diff suppressed because it is too large Load Diff

734
phono3py/conductivity/init_direct_solution.py Normal file
View File

@ -0,0 +1,734 @@
"""Init lattice thermal conductivity classes with direct solution."""
# Copyright (C) 2020 Atsushi Togo
# All rights reserved.
#
# This file is part of phono3py.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
#
# * Neither the name of the phonopy project nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import sys
from typing import Optional, Union
from phono3py.conductivity.base import get_unit_to_WmK
from phono3py.conductivity.direct_solution import ConductivityLBTE
from phono3py.conductivity.direct_solution_base import ConductivityLBTEBase
from phono3py.conductivity.utils import (
select_colmat_solver,
write_pp_interaction,
)
from phono3py.conductivity.wigner_direct_solution import ConductivityWignerLBTE
from phono3py.file_IO import (
read_collision_from_hdf5,
write_collision_eigenvalues_to_hdf5,
write_collision_to_hdf5,
write_kappa_to_hdf5,
write_unitary_matrix_to_hdf5,
)
from phono3py.phonon3.interaction import Interaction, all_bands_exist
cond_LBTE_type = Union[ConductivityLBTE, ConductivityWignerLBTE]
def get_thermal_conductivity_LBTE(
interaction: Interaction,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
is_isotope=False,
mass_variances=None,
grid_points=None,
boundary_mfp=None, # in micrometer
solve_collective_phonon=False,
is_reducible_collision_matrix=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
conductivity_type=None,
pinv_cutoff=1.0e-8,
pinv_solver=0, # default: dsyev in lapacke
pinv_method=0, # default: abs(eig) < cutoff
write_collision=False,
read_collision=False,
write_kappa=False,
write_pp=False,
read_pp=False,
write_LBTE_solution=False,
compression="gzip",
input_filename=None,
output_filename=None,
log_level=0,
):
"""Calculate lattice thermal conductivity by direct solution."""
if temperatures is None:
_temperatures = [
300,
]
else:
_temperatures = temperatures
if sigmas is None:
sigmas = []
if log_level:
print("-" * 19 + " Lattice thermal conductivity (LBTE) " + "-" * 19)
print(
"Cutoff frequency of pseudo inversion of collision matrix: %s" % pinv_cutoff
)
if read_collision:
temps = None
else:
temps = _temperatures
if conductivity_type == "wigner":
conductivity_LBTE_class = ConductivityWignerLBTE
else:
conductivity_LBTE_class = ConductivityLBTE
lbte = conductivity_LBTE_class(
interaction,
grid_points=grid_points,
temperatures=temps,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
boundary_mfp=boundary_mfp,
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=input_filename,
pinv_cutoff=pinv_cutoff,
pinv_solver=pinv_solver,
pinv_method=pinv_method,
log_level=log_level,
)
if read_collision:
read_from = _set_collision_from_file(
lbte,
indices=read_collision,
is_reducible_collision_matrix=is_reducible_collision_matrix,
filename=input_filename,
log_level=log_level,
)
if not read_from:
print("Reading collision failed.")
return False
if log_level:
temps_read = lbte.temperatures
if len(temps_read) > 5:
text = (" %.1f " * 5 + "...") % tuple(temps_read[:5])
text += " %.1f" % temps_read[-1]
else:
text = (" %.1f " * len(temps_read)) % tuple(temps_read)
print("Temperature: " + text)
# This computes pieces of collision matrix sequentially.
for i in lbte:
if write_pp:
write_pp_interaction(
lbte, interaction, i, filename=output_filename, compression=compression
)
if write_collision:
ConductivityLBTEWriter.write_collision(
lbte,
interaction,
i=i,
is_reducible_collision_matrix=is_reducible_collision_matrix,
is_one_gp_colmat=(grid_points is not None),
filename=output_filename,
)
lbte.delete_gp_collision_and_pp()
# Write full collision matrix
if write_LBTE_solution:
if (
read_collision
and all_bands_exist(interaction)
and read_from == "grid_points"
and grid_points is None
) or (not read_collision):
ConductivityLBTEWriter.write_collision(
lbte, interaction, filename=output_filename
)
if grid_points is None and all_bands_exist(interaction):
lbte.set_kappa_at_sigmas()
if write_kappa:
ConductivityLBTEWriter.write_kappa(
lbte,
interaction.primitive.volume,
is_reducible_collision_matrix=is_reducible_collision_matrix,
write_LBTE_solution=write_LBTE_solution,
pinv_solver=pinv_solver,
compression=compression,
filename=output_filename,
log_level=log_level,
)
return lbte
class ConductivityLBTEWriter:
"""Collection of result writers."""
@staticmethod
def write_collision(
lbte: cond_LBTE_type,
interaction: Interaction,
i=None,
is_reducible_collision_matrix=False,
is_one_gp_colmat=False,
filename=None,
):
"""Write collision matrix into hdf5 file."""
grid_points = lbte.grid_points
temperatures = lbte.temperatures
sigmas = lbte.sigmas
sigma_cutoff = lbte.sigma_cutoff_width
gamma = lbte.gamma
gamma_isotope = lbte.gamma_isotope
collision_matrix = lbte.collision_matrix
mesh = lbte.mesh_numbers
if i is not None:
gp = grid_points[i]
if is_one_gp_colmat:
igp = 0
else:
if is_reducible_collision_matrix:
igp = interaction.bz_grid.bzg2grg[gp]
else:
igp = i
if all_bands_exist(interaction):
for j, sigma in enumerate(sigmas):
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[j, igp]
else:
gamma_isotope_at_sigma = None
write_collision_to_hdf5(
temperatures,
mesh,
gamma=gamma[j, :, igp],
gamma_isotope=gamma_isotope_at_sigma,
collision_matrix=collision_matrix[j, :, igp],
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
)
else:
for j, sigma in enumerate(sigmas):
for k, bi in enumerate(interaction.band_indices):
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[j, igp, k]
else:
gamma_isotope_at_sigma = None
write_collision_to_hdf5(
temperatures,
mesh,
gamma=gamma[j, :, igp, k],
gamma_isotope=gamma_isotope_at_sigma,
collision_matrix=collision_matrix[j, :, igp, k],
grid_point=gp,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
)
else:
for j, sigma in enumerate(sigmas):
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[j]
else:
gamma_isotope_at_sigma = None
write_collision_to_hdf5(
temperatures,
mesh,
gamma=gamma[j],
gamma_isotope=gamma_isotope_at_sigma,
collision_matrix=collision_matrix[j],
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
)
@staticmethod
def write_kappa(
lbte: cond_LBTE_type,
volume: float,
is_reducible_collision_matrix: bool = False,
write_LBTE_solution: bool = False,
pinv_solver: Optional[int] = None,
compression: str = "gzip",
filename: Optional[str] = None,
log_level: int = 0,
):
"""Write kappa related properties into a hdf5 file."""
if isinstance(lbte, ConductivityLBTE):
kappa = lbte.kappa
mode_kappa = lbte.mode_kappa
kappa_RTA = lbte.kappa_RTA
mode_kappa_RTA = lbte.mode_kappa_RTA
gv = lbte.group_velocities
gv_by_gv = lbte.gv_by_gv
else:
kappa = None
mode_kappa = None
kappa_RTA = None
mode_kappa_RTA = None
gv = None
gv_by_gv = None
if isinstance(lbte, ConductivityWignerLBTE):
kappa_P_exact = lbte.kappa_P_exact
kappa_P_RTA = lbte.kappa_P_RTA
kappa_C = lbte.kappa_C
mode_kappa_P_exact = lbte.mode_kappa_P_exact
mode_kappa_P_RTA = lbte.mode_kappa_P_RTA
mode_kappa_C = lbte.mode_kappa_C
else:
kappa_P_exact = None
kappa_P_RTA = None
kappa_C = None
mode_kappa_P_exact = None
mode_kappa_P_RTA = None
mode_kappa_C = None
temperatures = lbte.temperatures
sigmas = lbte.sigmas
sigma_cutoff = lbte.sigma_cutoff_width
mesh = lbte.mesh_numbers
bz_grid = lbte.bz_grid
grid_points = lbte.grid_points
weights = lbte.grid_weights
frequencies = lbte.frequencies
ave_pp = lbte.averaged_pp_interaction
qpoints = lbte.qpoints
gamma = lbte.gamma
gamma_isotope = lbte.gamma_isotope
f_vector = lbte.get_f_vectors()
mode_cv = lbte.mode_heat_capacities
mfp = lbte.get_mean_free_path()
boundary_mfp = lbte.boundary_mfp
coleigs = lbte.collision_eigenvalues
# After kappa calculation, the variable is overwritten by unitary matrix
unitary_matrix = lbte.collision_matrix
if is_reducible_collision_matrix:
frequencies = lbte.get_frequencies_all()
else:
frequencies = lbte.frequencies
for i, sigma in enumerate(sigmas):
if kappa is None:
kappa_at_sigma = None
else:
kappa_at_sigma = kappa[i]
if mode_kappa is None:
mode_kappa_at_sigma = None
else:
mode_kappa_at_sigma = mode_kappa[i]
if kappa_RTA is None:
kappa_RTA_at_sigma = None
else:
kappa_RTA_at_sigma = kappa_RTA[i]
if mode_kappa_RTA is None:
mode_kappa_RTA_at_sigma = None
else:
mode_kappa_RTA_at_sigma = mode_kappa_RTA[i]
if kappa_P_exact is None:
kappa_P_exact_at_sigma = None
else:
kappa_P_exact_at_sigma = kappa_P_exact[i]
if kappa_P_RTA is None:
kappa_P_RTA_at_sigma = None
else:
kappa_P_RTA_at_sigma = kappa_P_RTA[i]
if kappa_C is None:
kappa_C_at_sigma = None
else:
kappa_C_at_sigma = kappa_C[i]
if kappa_P_exact_at_sigma is None or kappa_C_at_sigma is None:
kappa_TOT_exact_at_sigma = None
else:
kappa_TOT_exact_at_sigma = kappa_P_exact_at_sigma + kappa_C_at_sigma
if kappa_P_RTA_at_sigma is None or kappa_C_at_sigma is None:
kappa_TOT_RTA_at_sigma = None
else:
kappa_TOT_RTA_at_sigma = kappa_P_RTA_at_sigma + kappa_C_at_sigma
if mode_kappa_P_exact is None:
mode_kappa_P_exact_at_sigma = None
else:
mode_kappa_P_exact_at_sigma = mode_kappa_P_exact[i]
if mode_kappa_P_RTA is None:
mode_kappa_P_RTA_at_sigma = None
else:
mode_kappa_P_RTA_at_sigma = mode_kappa_P_RTA[i]
if mode_kappa_C is None:
mode_kappa_C_at_sigma = None
else:
mode_kappa_C_at_sigma = mode_kappa_C[i]
if gamma_isotope is None:
gamma_isotope_at_sigma = None
else:
gamma_isotope_at_sigma = gamma_isotope[i]
write_kappa_to_hdf5(
temperatures,
mesh,
boundary_mfp=boundary_mfp,
bz_grid=bz_grid,
frequency=frequencies,
group_velocity=gv,
gv_by_gv=gv_by_gv,
mean_free_path=mfp[i],
heat_capacity=mode_cv,
kappa=kappa_at_sigma,
mode_kappa=mode_kappa_at_sigma,
kappa_RTA=kappa_RTA_at_sigma,
mode_kappa_RTA=mode_kappa_RTA_at_sigma,
kappa_P_exact=kappa_P_exact_at_sigma,
kappa_P_RTA=kappa_P_RTA_at_sigma,
kappa_C=kappa_C_at_sigma,
kappa_TOT_exact=kappa_TOT_exact_at_sigma,
kappa_TOT_RTA=kappa_TOT_RTA_at_sigma,
mode_kappa_P_exact=mode_kappa_P_exact_at_sigma,
mode_kappa_P_RTA=mode_kappa_P_RTA_at_sigma,
mode_kappa_C=mode_kappa_C_at_sigma,
f_vector=f_vector,
gamma=gamma[i],
gamma_isotope=gamma_isotope_at_sigma,
averaged_pp_interaction=ave_pp,
qpoint=qpoints,
grid_point=grid_points,
weight=weights,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=get_unit_to_WmK() / volume,
compression=compression,
filename=filename,
verbose=log_level,
)
if coleigs is not None:
write_collision_eigenvalues_to_hdf5(
temperatures,
mesh,
coleigs[i],
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=log_level,
)
if write_LBTE_solution:
if pinv_solver is not None:
solver = select_colmat_solver(pinv_solver)
if solver in [1, 2, 3, 4, 5]:
write_unitary_matrix_to_hdf5(
temperatures,
mesh,
unitary_matrix=unitary_matrix,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
solver=solver,
filename=filename,
verbose=log_level,
)
def _set_collision_from_file(
lbte: ConductivityLBTEBase,
indices="all",
is_reducible_collision_matrix=False,
filename=None,
log_level=0,
):
"""Set collision matrix from that read from files.
If collision-m*.hdf5 that contains all data is not found,
collision-m*-gp*.hdf5 files at grid points are searched. If any of those
files are not found, collision-m*-gp*-b*.hdf5 files at grid points and bands
are searched. If any of those files are not found, it fails.
lbte : ConductivityLBTEBase
RTA lattice thermal conductivity instance.
filename : str, optional
This string is inserted in the filename as collision-m*.{filename}.hdf5.
verbose : bool, optional
Show text output or not.
"""
bz_grid = lbte.bz_grid
sigmas = lbte.sigmas
sigma_cutoff = lbte.sigma_cutoff_width
mesh = lbte.mesh_numbers
grid_points = lbte.grid_points
read_from = None
if log_level:
print(
"---------------------- Reading collision data from file "
"----------------------",
flush=True,
)
arrays_allocated = False
for i_sigma, sigma in enumerate(sigmas):
collisions = read_collision_from_hdf5(
mesh,
indices=indices,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=(log_level > 0),
)
if log_level:
sys.stdout.flush()
if collisions:
(colmat_at_sigma, gamma_at_sigma, temperatures) = collisions
if not arrays_allocated:
arrays_allocated = True
# The following invokes self._allocate_values()
lbte.temperatures = temperatures
lbte.collision_matrix[i_sigma] = colmat_at_sigma[0]
lbte.gamma[i_sigma] = gamma_at_sigma[0]
read_from = "full_matrix"
else:
vals = _allocate_collision(
True,
mesh,
sigma,
sigma_cutoff,
grid_points,
indices,
filename,
)
if vals:
colmat_at_sigma, gamma_at_sigma, temperatures = vals
else:
if log_level:
print("Collision at grid point %d doesn't exist." % grid_points[0])
vals = _allocate_collision(
False,
mesh,
sigma,
sigma_cutoff,
grid_points,
indices,
filename,
)
if vals:
colmat_at_sigma, gamma_at_sigma, temperatures = vals
else:
if log_level:
print(
"Collision at (grid point %d, band index %d) "
"doesn't exist." % (grid_points[0], 1)
)
return False
if not arrays_allocated:
arrays_allocated = True
# The following invokes self._allocate_values()
lbte.temperatures = temperatures
for i, gp in enumerate(grid_points):
if not _collect_collision_gp(
lbte.collision_matrix[i_sigma],
lbte.gamma[i_sigma],
temperatures,
mesh,
sigma,
sigma_cutoff,
i,
gp,
bz_grid.bzg2grg,
indices,
is_reducible_collision_matrix,
filename,
log_level,
):
num_band = lbte.collision_matrix.shape[3]
for i_band in range(num_band):
if not _collect_collision_band(
lbte.collision_matrix[i_sigma],
lbte.gamma[i_sigma],
temperatures,
mesh,
sigma,
sigma_cutoff,
i,
gp,
bz_grid.bzg2grg,
i_band,
indices,
is_reducible_collision_matrix,
filename,
log_level,
):
return False
read_from = "grid_points"
return read_from
def _allocate_collision(
for_gps,
mesh,
sigma,
sigma_cutoff,
grid_points,
indices,
filename,
):
if for_gps:
collision = read_collision_from_hdf5(
mesh,
indices=indices,
grid_point=grid_points[0],
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
only_temperatures=True,
verbose=False,
)
else:
collision = read_collision_from_hdf5(
mesh,
indices=indices,
grid_point=grid_points[0],
band_index=0,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
only_temperatures=True,
verbose=False,
)
if collision is None:
return False
temperatures = collision[2]
return None, None, temperatures
def _collect_collision_gp(
colmat_at_sigma,
gamma_at_sigma,
temperatures,
mesh,
sigma,
sigma_cutoff,
i,
gp,
bzg2grg,
indices,
is_reducible_collision_matrix,
filename,
log_level,
):
collision_gp = read_collision_from_hdf5(
mesh,
indices=indices,
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=(log_level > 0),
)
if log_level:
sys.stdout.flush()
if not collision_gp:
return False
(colmat_at_gp, gamma_at_gp, temperatures_at_gp) = collision_gp
if is_reducible_collision_matrix:
igp = bzg2grg[gp]
else:
igp = i
gamma_at_sigma[:, igp] = gamma_at_gp
colmat_at_sigma[:, igp] = colmat_at_gp[0]
temperatures[:] = temperatures_at_gp
return True
def _collect_collision_band(
colmat_at_sigma,
gamma_at_sigma,
temperatures,
mesh,
sigma,
sigma_cutoff,
i,
gp,
bzg2grg,
j,
indices,
is_reducible_collision_matrix,
filename,
log_level,
):
collision_band = read_collision_from_hdf5(
mesh,
indices=indices,
grid_point=gp,
band_index=j,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=(log_level > 0),
)
if log_level:
sys.stdout.flush()
if collision_band is False:
return False
(colmat_at_band, gamma_at_band, temperatures_at_band) = collision_band
if is_reducible_collision_matrix:
igp = bzg2grg[gp]
else:
igp = i
gamma_at_sigma[:, igp, j] = gamma_at_band[0]
colmat_at_sigma[:, igp, j] = colmat_at_band[0]
temperatures[:] = temperatures_at_band
return True

735
phono3py/conductivity/init_rta.py Normal file
View File

@ -0,0 +1,735 @@
"""Init lattice thermal conductivity classes with RTA."""
# Copyright (C) 2020 Atsushi Togo
# All rights reserved.
#
# This file is part of phono3py.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
#
# * Neither the name of the phonopy project nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from typing import Optional, Union, cast
import numpy as np
from phono3py.conductivity.base import get_unit_to_WmK
from phono3py.conductivity.kubo_rta import ConductivityKuboRTA
from phono3py.conductivity.rta import ConductivityRTA
from phono3py.conductivity.rta_base import ConductivityRTABase
from phono3py.conductivity.utils import write_pp_interaction
from phono3py.conductivity.wigner_rta import ConductivityWignerRTA
from phono3py.file_IO import (
read_gamma_from_hdf5,
write_gamma_detail_to_hdf5,
write_kappa_to_hdf5,
)
from phono3py.phonon3.interaction import Interaction, all_bands_exist
from phono3py.phonon3.triplets import get_all_triplets
cond_RTA_type = Union[ConductivityRTA, ConductivityWignerRTA, ConductivityKuboRTA]
def get_thermal_conductivity_RTA(
interaction: Interaction,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
mass_variances=None,
grid_points=None,
is_isotope=False,
boundary_mfp=None, # in micrometer
use_ave_pp=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
is_N_U=False,
conductivity_type=None,
write_gamma=False,
read_gamma=False,
write_kappa=False,
write_pp=False,
read_pp=False,
write_gamma_detail=False,
compression="gzip",
input_filename=None,
output_filename=None,
log_level=0,
):
"""Run RTA thermal conductivity calculation."""
if temperatures is None:
_temperatures = np.arange(0, 1001, 10, dtype="double")
else:
_temperatures = temperatures
if conductivity_type == "wigner":
conductivity_RTA_class = ConductivityWignerRTA
elif conductivity_type == "kubo":
conductivity_RTA_class = ConductivityKuboRTA
else:
conductivity_RTA_class = ConductivityRTA
if log_level:
print(
"-------------------- Lattice thermal conductivity (RTA) "
"--------------------"
)
br = conductivity_RTA_class(
interaction,
grid_points=grid_points,
temperatures=_temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
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=write_pp,
pp_filename=input_filename,
is_N_U=is_N_U,
is_gamma_detail=write_gamma_detail,
log_level=log_level,
)
if read_gamma:
if not _set_gamma_from_file(br, filename=input_filename):
print("Reading collisions failed.")
return False
for i in br:
if write_pp:
write_pp_interaction(
br, interaction, i, compression=compression, filename=output_filename
)
if write_gamma:
ConductivityRTAWriter.write_gamma(
br,
interaction,
i,
compression=compression,
filename=output_filename,
verbose=log_level,
)
if write_gamma_detail:
ConductivityRTAWriter.write_gamma_detail(
br,
interaction,
i,
compression=compression,
filename=output_filename,
verbose=log_level,
)
if grid_points is None and all_bands_exist(interaction):
br.set_kappa_at_sigmas()
if log_level:
if conductivity_type == "wigner":
ShowCalcProgress.kappa_Wigner_RTA(
cast(ConductivityWignerRTA, br), log_level
)
else:
ShowCalcProgress.kappa_RTA(cast(ConductivityRTA, br), log_level)
if write_kappa:
ConductivityRTAWriter.write_kappa(
br,
interaction.primitive.volume,
compression=compression,
filename=output_filename,
log_level=log_level,
)
return br
class ShowCalcProgress:
"""Show calculation progress."""
@staticmethod
def kappa_RTA(br: ConductivityRTA, log_level):
"""Show RTA calculation progress."""
temperatures = br.temperatures
sigmas = br.sigmas
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
for i, sigma in enumerate(sigmas):
text = "----------- Thermal conductivity (W/m-k) "
if sigma:
text += "for sigma=%s -----------" % sigma
else:
text += "with tetrahedron method -----------"
print(text)
if log_level > 1:
print(
("#%6s " + " %-10s" * 6 + "#ipm")
% ("T(K)", "xx", "yy", "zz", "yz", "xz", "xy")
)
for j, (t, k) in enumerate(zip(temperatures, kappa[i])):
print(
("%7.1f" + " %10.3f" * 6 + " %d/%d")
% (
(t,)
+ tuple(k)
+ (num_ignored_phonon_modes[i, j], num_phonon_modes)
)
)
else:
print(
("#%6s " + " %-10s" * 6)
% ("T(K)", "xx", "yy", "zz", "yz", "xz", "xy")
)
for t, k in zip(temperatures, kappa[i]):
print(("%7.1f " + " %10.3f" * 6) % ((t,) + tuple(k)))
print("", flush=True)
@staticmethod
def kappa_Wigner_RTA(br: ConductivityWignerRTA, log_level):
"""Show Wigner-RTA calculation progress."""
temperatures = br.temperatures
sigmas = br.sigmas
kappa_TOT_RTA = br.kappa_TOT_RTA
kappa_P_RTA = br.kappa_P_RTA
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
for i, sigma in enumerate(sigmas):
text = "----------- Thermal conductivity (W/m-k) "
if sigma:
text += "for sigma=%s -----------" % sigma
else:
text += "with tetrahedron method -----------"
print(text)
if log_level > 1:
print(
("#%6s " + " %-10s" * 6 + "#ipm")
% (" \t T(K)", "xx", "yy", "zz", "yz", "xz", "xy")
)
for j, (t, k) in enumerate(zip(temperatures, kappa_P_RTA[i])):
print(
"K_P\t"
+ ("%7.1f" + " %10.3f" * 6 + " %d/%d")
% (
(t,)
+ tuple(k)
+ (num_ignored_phonon_modes[i, j], num_phonon_modes)
)
)
print(" ")
for j, (t, k) in enumerate(zip(temperatures, kappa_C[i])):
print(
"K_C\t"
+ ("%7.1f" + " %10.3f" * 6 + " %d/%d")
% (
(t,)
+ tuple(k)
+ (num_ignored_phonon_modes[i, j], num_phonon_modes)
)
)
print(" ")
for j, (t, k) in enumerate(zip(temperatures, kappa_TOT_RTA[i])):
print(
"K_T\t"
+ ("%7.1f" + " %10.3f" * 6 + " %d/%d")
% (
(t,)
+ tuple(k)
+ (num_ignored_phonon_modes[i, j], num_phonon_modes)
)
)
else:
print(
("#%6s " + " %-10s" * 6)
% (" \t T(K)", "xx", "yy", "zz", "yz", "xz", "xy")
)
if kappa_P_RTA is not None:
for t, k in zip(temperatures, kappa_P_RTA[i]):
print("K_P\t" + ("%7.1f " + " %10.3f" * 6) % ((t,) + tuple(k)))
print(" ")
for t, k in zip(temperatures, kappa_C[i]):
print("K_C\t" + ("%7.1f " + " %10.3f" * 6) % ((t,) + tuple(k)))
print(" ")
for t, k in zip(temperatures, kappa_TOT_RTA[i]):
print("K_T\t" + ("%7.1f " + " %10.3f" * 6) % ((t,) + tuple(k)))
print("", flush=True)
class ConductivityRTAWriter:
"""Collection of result writers."""
@staticmethod
def write_gamma(
br: cond_RTA_type,
interaction: Interaction,
i: int,
compression: str = "gzip",
filename: Optional[str] = None,
verbose: bool = True,
):
"""Write mode kappa related properties into a hdf5 file."""
grid_points = br.grid_points
if isinstance(br, ConductivityRTA):
group_velocities_i = br.group_velocities[i]
gv_by_gv_i = br.gv_by_gv[i]
else:
group_velocities_i = None
gv_by_gv_i = None
if isinstance(br, ConductivityWignerRTA):
velocity_operator_i = br.velocity_operator[i]
else:
velocity_operator_i = None
if isinstance(br, (ConductivityRTA, ConductivityWignerRTA)):
mode_heat_capacities = br.mode_heat_capacities
else:
mode_heat_capacities = None
ave_pp = br.averaged_pp_interaction
mesh = br.mesh_numbers
bz_grid = br.bz_grid
temperatures = br.temperatures
gamma = br.gamma
gamma_isotope = br.gamma_isotope
sigmas = br.sigmas
sigma_cutoff = br.sigma_cutoff_width
volume = interaction.primitive.volume
gamma_N, gamma_U = br.get_gamma_N_U()
gp = grid_points[i]
if all_bands_exist(interaction):
if ave_pp is None:
ave_pp_i = None
else:
ave_pp_i = ave_pp[i]
frequencies = interaction.get_phonons()[0][gp]
for j, sigma in enumerate(sigmas):
if gamma_isotope is None:
gamma_isotope_at_sigma = None
else:
gamma_isotope_at_sigma = gamma_isotope[j, i]
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[j, :, i]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[j, :, i]
write_kappa_to_hdf5(
temperatures,
mesh,
bz_grid=bz_grid,
frequency=frequencies,
group_velocity=group_velocities_i,
gv_by_gv=gv_by_gv_i,
velocity_operator=velocity_operator_i,
heat_capacity=mode_heat_capacities[:, i],
gamma=gamma[j, :, i],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp_i,
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=get_unit_to_WmK() / volume,
compression=compression,
filename=filename,
verbose=verbose,
)
else:
for j, sigma in enumerate(sigmas):
for k, bi in enumerate(interaction.band_indices):
if group_velocities_i is None:
group_velocities_ik = None
else:
group_velocities_ik = group_velocities_i[k]
if velocity_operator_i is None:
velocity_operator_ik = None
else:
velocity_operator_ik = velocity_operator_i[k]
if gv_by_gv_i is None:
gv_by_gv_ik = None
else:
gv_by_gv_ik = gv_by_gv_i[k]
if ave_pp is None:
ave_pp_ik = None
else:
ave_pp_ik = ave_pp[i, k]
frequencies = interaction.get_phonons()[0][gp, bi]
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[j, i, k]
else:
gamma_isotope_at_sigma = None
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[j, :, i, k]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[j, :, i, k]
write_kappa_to_hdf5(
temperatures,
mesh,
bz_grid=bz_grid,
frequency=frequencies,
group_velocity=group_velocities_ik,
gv_by_gv=gv_by_gv_ik,
velocity_operator=velocity_operator_ik,
heat_capacity=mode_heat_capacities[:, i, k],
gamma=gamma[j, :, i, k],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp_ik,
grid_point=gp,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=get_unit_to_WmK() / volume,
compression=compression,
filename=filename,
verbose=verbose,
)
@staticmethod
def write_kappa(
br: cond_RTA_type,
volume: float,
compression: str = "gzip",
filename: Optional[str] = None,
log_level: int = 0,
):
"""Write kappa related properties into a hdf5 file."""
temperatures = br.temperatures
sigmas = br.sigmas
sigma_cutoff = br.sigma_cutoff_width
gamma = br.gamma
gamma_isotope = br.gamma_isotope
gamma_N, gamma_U = br.get_gamma_N_U()
mesh = br.mesh_numbers
bz_grid = br.bz_grid
frequencies = br.frequencies
if isinstance(br, ConductivityRTA):
kappa = br.kappa
mode_kappa = br.mode_kappa
gv = br.group_velocities
gv_by_gv = br.gv_by_gv
else:
kappa = None
mode_kappa = None
gv = None
gv_by_gv = None
if isinstance(br, ConductivityWignerRTA):
kappa_TOT_RTA = br.kappa_TOT_RTA
kappa_P_RTA = br.kappa_P_RTA
kappa_C = br.kappa_C
mode_kappa_P_RTA = br.mode_kappa_P_RTA
mode_kappa_C = br.mode_kappa_C
else:
kappa_TOT_RTA = None
kappa_P_RTA = None
kappa_C = None
mode_kappa_P_RTA = None
mode_kappa_C = None
if isinstance(br, (ConductivityRTA, ConductivityWignerRTA)):
mode_cv = br.mode_heat_capacities
else:
mode_cv = None
ave_pp = br.averaged_pp_interaction
qpoints = br.qpoints
grid_points = br.grid_points
weights = br.grid_weights
boundary_mfp = br.boundary_mfp
for i, sigma in enumerate(sigmas):
if kappa is None:
kappa_at_sigma = None
else:
kappa_at_sigma = kappa[i]
if mode_kappa is None:
mode_kappa_at_sigma = None
else:
mode_kappa_at_sigma = mode_kappa[i]
if kappa_TOT_RTA is None:
kappa_TOT_RTA_at_sigma = None
else:
kappa_TOT_RTA_at_sigma = kappa_TOT_RTA[i]
if kappa_P_RTA is None:
kappa_P_RTA_at_sigma = None
else:
kappa_P_RTA_at_sigma = kappa_P_RTA[i]
if kappa_C is None:
kappa_C_at_sigma = None
else:
kappa_C_at_sigma = kappa_C[i]
if mode_kappa_P_RTA is None:
mode_kappa_P_RTA_at_sigma = None
else:
mode_kappa_P_RTA_at_sigma = mode_kappa_P_RTA[i]
if mode_kappa_C is None:
mode_kappa_C_at_sigma = None
else:
mode_kappa_C_at_sigma = mode_kappa_C[i]
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[i]
else:
gamma_isotope_at_sigma = None
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[i]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[i]
write_kappa_to_hdf5(
temperatures,
mesh,
boundary_mfp=boundary_mfp,
bz_grid=bz_grid,
frequency=frequencies,
group_velocity=gv,
gv_by_gv=gv_by_gv,
heat_capacity=mode_cv,
kappa=kappa_at_sigma,
mode_kappa=mode_kappa_at_sigma,
kappa_TOT_RTA=kappa_TOT_RTA_at_sigma,
kappa_P_RTA=kappa_P_RTA_at_sigma,
kappa_C=kappa_C_at_sigma,
mode_kappa_P_RTA=mode_kappa_P_RTA_at_sigma,
mode_kappa_C=mode_kappa_C_at_sigma,
gamma=gamma[i],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp,
qpoint=qpoints,
grid_point=grid_points,
weight=weights,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=get_unit_to_WmK() / volume,
compression=compression,
filename=filename,
verbose=log_level,
)
@staticmethod
def write_gamma_detail(
br: cond_RTA_type,
interaction: Interaction,
i: int,
compression: str = "gzip",
filename: Optional[str] = None,
verbose: bool = True,
):
"""Write detailed Gamma values to hdf5 files."""
gamma_detail = br.get_gamma_detail_at_q()
temperatures = br.temperatures
mesh = br.mesh_numbers
bz_grid = br.bz_grid
grid_points = br.grid_points
gp = grid_points[i]
sigmas = br.sigmas
sigma_cutoff = br.sigma_cutoff_width
triplets, weights, _, _ = interaction.get_triplets_at_q()
all_triplets = get_all_triplets(gp, interaction.bz_grid)
if all_bands_exist(interaction):
for sigma in sigmas:
write_gamma_detail_to_hdf5(
temperatures,
mesh,
bz_grid=bz_grid,
gamma_detail=gamma_detail,
grid_point=gp,
triplet=triplets,
weight=weights,
triplet_all=all_triplets,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
compression=compression,
filename=filename,
verbose=verbose,
)
else:
for sigma in sigmas:
for k, bi in enumerate(interaction.get_band_indices()):
write_gamma_detail_to_hdf5(
temperatures,
mesh,
bz_grid=bz_grid,
gamma_detail=gamma_detail[:, :, k, :, :],
grid_point=gp,
triplet=triplets,
weight=weights,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
compression=compression,
filename=filename,
verbose=verbose,
)
def _set_gamma_from_file(
br: ConductivityRTABase, filename: Optional[str] = None, verbose: bool = True
):
"""Read kappa-*.hdf5 files for thermal conductivity calculation.
If kappa-m*.hdf5 that contains all data is not found, kappa-m*-gp*.hdf5
files at grid points are searched. If any of those files are not found,
kappa-m*-gp*-b*.hdf5 files at grid points and bands are searched. If any
of those files are not found, it fails.
br : ConductivityRTABase
RTA lattice thermal conductivity instance.
filename : str, optional
This string is inserted in the filename as kappa-m*.{filename}.hdf5.
verbose : bool, optional
Show text output or not.
"""
sigmas = br.sigmas
sigma_cutoff = br.sigma_cutoff_width
mesh = br.mesh_numbers
grid_points = br.grid_points
temperatures = br.temperatures
num_band = br.frequencies.shape[1]
gamma = np.zeros(
(len(sigmas), len(temperatures), len(grid_points), num_band), dtype="double"
)
gamma_N = np.zeros_like(gamma)
gamma_U = np.zeros_like(gamma)
gamma_iso = np.zeros((len(sigmas), len(grid_points), num_band), dtype="double")
ave_pp = np.zeros((len(grid_points), num_band), dtype="double")
is_gamma_N_U_in = False
is_ave_pp_in = False
read_succeeded = True
for j, sigma in enumerate(sigmas):
data, full_filename = read_gamma_from_hdf5(
mesh,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
)
if data:
if verbose:
print("Read data from %s." % full_filename)
gamma[j] = data["gamma"]
if "gamma_isotope" in data:
gamma_iso[j] = data["gamma_isotope"]
if "gamma_N" in data:
is_gamma_N_U_in = True
gamma_N[j] = data["gamma_N"]
gamma_U[j] = data["gamma_U"]
if "ave_pp" in data:
is_ave_pp_in = True
ave_pp[:] = data["ave_pp"]
else:
if verbose:
print(
"%s not found. Look for hdf5 files at grid points." % full_filename
)
for i, gp in enumerate(grid_points):
data_gp, full_filename = read_gamma_from_hdf5(
mesh,
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
)
if data_gp:
if verbose:
print("Read data from %s." % full_filename)
gamma[j, :, i] = data_gp["gamma"]
if "gamma_iso" in data_gp:
gamma_iso[j, i] = data_gp["gamma_iso"]
if "gamma_N" in data_gp:
is_gamma_N_U_in = True
gamma_N[j, :, i] = data_gp["gamma_N"]
gamma_U[j, :, i] = data_gp["gamma_U"]
if "ave_pp" in data_gp:
is_ave_pp_in = True
ave_pp[i] = data_gp["ave_pp"]
else:
if verbose:
print(
"%s not found. Look for hdf5 files at bands."
% full_filename
)
for bi in range(num_band):
data_band, full_filename = read_gamma_from_hdf5(
mesh,
grid_point=gp,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
)
if data_band:
if verbose:
print("Read data from %s." % full_filename)
gamma[j, :, i, bi] = data_band["gamma"]
if "gamma_iso" in data_band:
gamma_iso[j, i, bi] = data_band["gamma_iso"]
if "gamma_N" in data_band:
is_gamma_N_U_in = True
gamma_N[j, :, i, bi] = data_band["gamma_N"]
gamma_U[j, :, i, bi] = data_band["gamma_U"]
if "ave_pp" in data_band:
is_ave_pp_in = True
ave_pp[i, bi] = data_band["ave_pp"]
else:
if verbose:
print("%s not found." % full_filename)
read_succeeded = False
if read_succeeded:
br.gamma = gamma
if is_ave_pp_in:
br.set_averaged_pp_interaction(ave_pp)
if is_gamma_N_U_in:
br.set_gamma_N_U(gamma_N, gamma_U)
return True
else:
return False

8
phono3py/conductivity/kubo.py → phono3py/conductivity/kubo_base.py
View File

@ -35,7 +35,7 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.units import Kb, THzToEv
from phonopy.physical_units import get_physical_units
from phono3py.phonon.group_velocity_matrix import GroupVelocityMatrix
from phono3py.phonon.heat_capacity_matrix import mode_cv_matrix
@ -71,11 +71,11 @@ class ConductivityKuboMixIn:
"""
irgp = self._grid_points[i_gp]
freqs = self._frequencies[irgp] * THzToEv
cutoff = self._pp.cutoff_frequency * THzToEv
freqs = self._frequencies[irgp] * get_physical_units().THzToEv
cutoff = self._pp.cutoff_frequency * get_physical_units().THzToEv
for i_temp, temp in enumerate(self._temperatures):
if (freqs / (temp * Kb) > 100).any():
if (freqs / (temp * get_physical_units().KB) > 100).any():
continue
cvm = mode_cv_matrix(temp, freqs, cutoff=cutoff)
self._cv_mat[i_temp, i_data] = cvm[self._pp.band_indices, :]

188
phono3py/conductivity/kubo_rta.py Normal file
View File

@ -0,0 +1,188 @@
"""Kubo thermal conductivity RTA class."""
# Copyright (C) 2022 Atsushi Togo
# All rights reserved.
#
# This file is part of phono3py.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
#
# * Neither the name of the phonopy project nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phono3py.conductivity.kubo_base import ConductivityKuboMixIn
from phono3py.conductivity.rta_base import ConductivityRTABase
from phono3py.phonon3.interaction import Interaction
class ConductivityKuboRTA(ConductivityKuboMixIn, ConductivityRTABase):
"""Class of Kubo lattice thermal conductivity under 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,
):
"""Init method."""
self._cv_mat = None
self._gv_mat = None
self._kappa = None
super().__init__(
interaction,
grid_points=grid_points,
temperatures=temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
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,
pp_filename=pp_filename,
is_N_U=is_N_U,
is_gamma_detail=is_gamma_detail,
is_frequency_shift_by_bubble=is_frequency_shift_by_bubble,
log_level=log_level,
)
def set_kappa_at_sigmas(self):
"""Calculate kappa from ph-ph interaction results."""
for i_gp, _ in enumerate(self._grid_points):
frequencies = self._frequencies[self._grid_points[i_gp]]
for j in range(len(self._sigmas)):
for k in range(len(self._temperatures)):
g_sum = self._get_main_diagonal(i_gp, j, k)
for i_band, freq in enumerate(frequencies):
if freq < self._pp.cutoff_frequency:
self._num_ignored_phonon_modes[j, k] += 1
continue
self._set_kappa_at_sigmas(
j, k, i_gp, i_band, g_sum, frequencies
)
N = self._num_sampling_grid_points
self._kappa = self._mode_kappa_mat.sum(axis=2).sum(axis=2).sum(axis=2).real / N
def _set_kappa_at_sigmas(self, j, k, i_gp, i_band, g_sum, frequencies):
gvm_sum2 = self._gv_mat_sum2[i_gp]
cvm = self._cv_mat[k, i_gp]
for j_band, freq in enumerate(frequencies):
if freq < self._pp.cutoff_frequency:
return
g = g_sum[i_band] + g_sum[j_band]
for i_pair, _ in enumerate(
([0, 0], [1, 1], [2, 2], [1, 2], [0, 2], [0, 1])
):
old_settings = np.seterr(all="raise")
try:
self._mode_kappa_mat[j, k, i_gp, i_band, j_band, i_pair] = (
cvm[i_band, j_band]
* gvm_sum2[i_band, j_band, i_pair]
* g
/ ((frequencies[j_band] - frequencies[i_band]) ** 2 + g**2)
* self._conversion_factor
)
except FloatingPointError:
# supposed that g is almost 0 and |gv|=0
pass
except Exception:
gp = self._grid_points[i_gp]
print("=" * 26 + " Warning " + "=" * 26)
print(
" Unexpected physical condition of ph-ph "
"interaction calculation was found."
)
print(
" g=%f at gp=%d, band=%d, freq=%f, band=%d, freq=%f"
% (
g_sum[i_band],
gp,
i_band + 1,
frequencies[i_band],
j_band + 1,
frequencies[j_band],
)
)
print("=" * 61)
np.seterr(**old_settings)
def _allocate_values(self):
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_mat = np.zeros(
(num_temp, num_grid_points, num_band0, num_band), dtype="double", order="C"
)
self._gv_mat = np.zeros(
(num_grid_points, num_band0, num_band, 3),
dtype=self._complex_dtype,
order="C",
)
self._gv_mat_sum2 = np.zeros(
(num_grid_points, num_band0, num_band, 6),
dtype=self._complex_dtype,
order="C",
)
# kappa and mode_kappa_mat are accessed when all bands exist, i.e.,
# num_band0==num_band.
self._kappa = np.zeros(
(len(self._sigmas), num_temp, 6), dtype="double", order="C"
)
self._mode_kappa_mat = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band, num_band, 6),
dtype=self._complex_dtype,
order="C",
)

962
phono3py/conductivity/rta.py
View File

@ -34,497 +34,11 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import warnings
from abc import abstractmethod
from typing import Type, Union, cast
import numpy as np
from phonopy.units import THzToEv
from phono3py.conductivity.base import ConductivityBase, ConductivityMixIn
from phono3py.conductivity.kubo import ConductivityKuboMixIn
from phono3py.conductivity.utils import (
ConductivityRTAWriter,
ShowCalcProgress,
set_gamma_from_file,
)
from phono3py.conductivity.utils import write_pp as write_phph
from phono3py.conductivity.wigner import (
ConductivityWignerMixIn,
get_conversion_factor_WTE,
)
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.phonon3.imag_self_energy import ImagSelfEnergy, average_by_degeneracy
from phono3py.phonon3.interaction import Interaction, all_bands_exist
class ConductivityRTABase(ConductivityBase):
"""Base class of ConductivityRTA*.
This is a base class of RTA classes.
"""
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,
):
"""Init method."""
self._is_N_U = is_N_U
self._is_gamma_detail = is_gamma_detail
self._is_frequency_shift_by_bubble = is_frequency_shift_by_bubble
self._gamma_N = None
self._gamma_U = None
self._gamma_detail_at_q = None
self._use_ave_pp = use_ave_pp
self._use_const_ave_pp = None
self._num_ignored_phonon_modes = None
super().__init__(
interaction,
grid_points=grid_points,
temperatures=temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
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,
)
self._use_const_ave_pp = self._pp.constant_averaged_interaction
self._read_pp = read_pp
self._store_pp = store_pp
self._pp_filename = pp_filename
if self._temperatures is not None:
self._allocate_values()
self._collision = ImagSelfEnergy(
self._pp, with_detail=(self._is_gamma_detail or self._is_N_U)
)
def get_gamma_N_U(self):
"""Return N and U parts of gamma."""
return (self._gamma_N, self._gamma_U)
def set_gamma_N_U(self, gamma_N, gamma_U):
"""Set N and U parts of gamma."""
self._gamma_N = gamma_N
self._gamma_U = gamma_U
def get_gamma_detail_at_q(self):
"""Return contribution of each triplet to gamma at current q-point."""
return self._gamma_detail_at_q
def get_number_of_ignored_phonon_modes(self):
"""Return number of ignored phonon modes."""
warnings.warn(
"Use attribute, number_of_ignored_phonon_modes",
DeprecationWarning,
stacklevel=2,
)
return self.number_of_ignored_phonon_modes
@property
def number_of_ignored_phonon_modes(self):
"""Return number of ignored phonon modes."""
return self._num_ignored_phonon_modes
def set_averaged_pp_interaction(self, ave_pp):
"""Set averaged ph-ph interaction."""
self._averaged_pp_interaction = ave_pp
@abstractmethod
def set_kappa_at_sigmas(self):
"""Must be implemented in the inherited class."""
raise NotImplementedError()
def _allocate_values(self):
num_band0 = len(self._pp.band_indices)
num_grid_points = len(self._grid_points)
num_temp = len(self._temperatures)
if not self._read_gamma:
self._gamma = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band0),
order="C",
dtype="double",
)
if self._is_gamma_detail or self._is_N_U:
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),
order="C",
dtype="double",
)
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
self._averaged_pp_interaction = np.zeros(
(num_grid_points, num_band0), order="C", dtype="double"
)
self._num_ignored_phonon_modes = np.zeros(
(len(self._sigmas), num_temp), order="C", dtype="intc"
)
def _run_at_grid_point(self):
i_gp = self._grid_point_count
self._show_log_header(i_gp)
grid_point = self._grid_points[i_gp]
self._set_cv(i_gp, i_gp)
self._set_velocities(i_gp, i_gp)
if self._read_gamma:
if self._use_ave_pp:
self._collision.set_grid_point(grid_point)
self._set_gamma_at_sigmas(i_gp)
else:
self._collision.set_grid_point(grid_point)
num_triplets = len(self._pp.get_triplets_at_q()[0])
if self._log_level:
print("Number of triplets: %d" % num_triplets)
if (
self._is_full_pp
or self._read_pp
or self._store_pp
or self._use_ave_pp
or self._use_const_ave_pp
or self._is_gamma_detail
):
self._set_gamma_at_sigmas(i_gp)
else: # can save memory space
self._set_gamma_at_sigmas_lowmem(i_gp)
if self._is_isotope and not self._read_gamma_iso:
gamma_iso = self._get_gamma_isotope_at_sigmas(i_gp)
self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._pp.band_indices]
if self._log_level:
self._show_log(i_gp)
def _set_gamma_at_sigmas(self, i):
for j, sigma in enumerate(self._sigmas):
self._collision.set_sigma(sigma, sigma_cutoff=self._sigma_cutoff)
self._collision.run_integration_weights()
if self._log_level:
text = "Collisions will be calculated with "
if sigma is None:
text += "tetrahedron method."
else:
text += "sigma=%s" % sigma
if self._sigma_cutoff is None:
text += "."
else:
text += "(%4.2f SD)." % self._sigma_cutoff
print(text)
if self._read_pp:
pp, _g_zero = read_pp_from_hdf5(
self._pp.mesh_numbers,
grid_point=self._grid_points[i],
sigma=sigma,
sigma_cutoff=self._sigma_cutoff,
filename=self._pp_filename,
verbose=(self._log_level > 0),
)
_, g_zero = self._collision.get_integration_weights()
if self._log_level:
if len(self._sigmas) > 1:
print(
"Multiple sigmas or mixing smearing and "
"tetrahedron method is not supported."
)
if _g_zero is not None and (_g_zero != g_zero).any():
raise ValueError("Inconsistency found in g_zero.")
self._collision.set_interaction_strength(pp)
elif self._use_ave_pp:
self._collision.set_averaged_pp_interaction(
self._averaged_pp_interaction[i]
)
elif self._use_const_ave_pp:
if self._log_level:
print(
"Constant ph-ph interaction of %6.3e is used."
% self._pp.constant_averaged_interaction
)
self._collision.run_interaction()
self._averaged_pp_interaction[i] = self._pp.averaged_interaction
elif j != 0 and (self._is_full_pp or self._sigma_cutoff is None):
if self._log_level:
print("Existing ph-ph interaction is used.")
else:
if self._log_level:
print("Calculating ph-ph interaction...")
self._collision.run_interaction(is_full_pp=self._is_full_pp)
if self._is_full_pp:
self._averaged_pp_interaction[i] = self._pp.averaged_interaction
# Number of triplets depends on q-point.
# So this is allocated each time.
if self._is_gamma_detail:
num_temp = len(self._temperatures)
self._gamma_detail_at_q = np.empty(
((num_temp,) + self._pp.get_interaction_strength().shape),
dtype="double",
order="C",
)
self._gamma_detail_at_q[:] = 0
if self._log_level:
print("Calculating collisions at temperatures...")
for k, t in enumerate(self._temperatures):
self._collision.temperature = t
self._collision.run()
self._gamma[j, k, i] = self._collision.imag_self_energy
if self._is_N_U or self._is_gamma_detail:
g_N, g_U = self._collision.get_imag_self_energy_N_and_U()
self._gamma_N[j, k, i] = g_N
self._gamma_U[j, k, i] = g_U
if self._is_gamma_detail:
self._gamma_detail_at_q[k] = (
self._collision.get_detailed_imag_self_energy()
)
def _set_gamma_at_sigmas_lowmem(self, i):
"""Calculate gamma without storing ph-ph interaction strength.
`svecs` and `multi` below must not be simply replaced by
`self._pp.primitive.get_smallest_vectors()` because they must be in
dense format as always so in Interaction class instance.
`p2s`, `s2p`, and `masses` have to be also given from Interaction
class instance.
"""
num_band = len(self._pp.primitive) * 3
band_indices = self._pp.band_indices
(
svecs,
multi,
p2s,
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:
tetrahedra = get_tetrahedra_relative_grid_address(
self._pp.bz_grid.microzone_lattice
)
# It is assumed that self._sigmas = [None].
for j, sigma in enumerate(self._sigmas):
self._collision.set_sigma(sigma)
if self._is_N_U:
collisions = np.zeros(
(2, len(self._temperatures), len(band_indices)),
dtype="double",
order="C",
)
else:
collisions = np.zeros(
(len(self._temperatures), len(band_indices)),
dtype="double",
order="C",
)
import phono3py._phono3py as phono3c
# True: OpenMP over triplets
# False: OpenMP over bands
if self._pp.openmp_per_triplets is None:
if len(triplets_at_q) > num_band:
openmp_per_triplets = True
else:
openmp_per_triplets = False
else:
openmp_per_triplets = self._pp.openmp_per_triplets
if sigma is None:
phono3c.pp_collision(
collisions,
np.array(
np.dot(tetrahedra, self._pp.bz_grid.P.T),
dtype="int64",
order="C",
),
self._frequencies,
self._eigenvectors,
triplets_at_q,
weights_at_q,
self._pp.bz_grid.addresses,
self._pp.bz_grid.gp_map,
self._pp.bz_grid.store_dense_gp_map * 1 + 1,
self._pp.bz_grid.D_diag,
self._pp.bz_grid.Q,
fc3,
svecs,
multi,
masses,
p2s,
s2p,
band_indices,
self._temperatures,
self._is_N_U * 1,
self._pp.symmetrize_fc3q * 1,
self._pp.make_r0_average * 1,
self._pp.all_shortest,
self._pp.cutoff_frequency,
openmp_per_triplets * 1,
)
else:
if self._sigma_cutoff is None:
sigma_cutoff = -1
else:
sigma_cutoff = float(self._sigma_cutoff)
phono3c.pp_collision_with_sigma(
collisions,
sigma,
sigma_cutoff,
self._frequencies,
self._eigenvectors,
triplets_at_q,
weights_at_q,
self._pp.bz_grid.addresses,
self._pp.bz_grid.D_diag,
self._pp.bz_grid.Q,
fc3,
svecs,
multi,
masses,
p2s,
s2p,
band_indices,
self._temperatures,
self._is_N_U * 1,
self._pp.symmetrize_fc3q * 1,
self._pp.make_r0_average * 1,
self._pp.all_shortest,
self._pp.cutoff_frequency,
openmp_per_triplets * 1,
)
col_unit_conv = self._collision.unit_conversion_factor
pp_unit_conv = self._pp.unit_conversion_factor
if self._is_N_U:
col = collisions.sum(axis=0)
col_N = collisions[0]
col_U = collisions[1]
else:
col = collisions
for k in range(len(self._temperatures)):
self._gamma[j, k, i, :] = average_by_degeneracy(
col[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]],
)
if self._is_N_U:
self._gamma_N[j, k, i, :] = average_by_degeneracy(
col_N[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]],
)
self._gamma_U[j, k, i, :] = average_by_degeneracy(
col_U[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]],
)
def _show_log(self, i_gp):
q = self._get_qpoint_from_gp_index(i_gp)
gp = self._grid_points[i_gp]
frequencies = self._frequencies[gp][self._pp.band_indices]
gv = self._gv[i_gp]
if self._averaged_pp_interaction is not None:
ave_pp = self._averaged_pp_interaction[i_gp]
else:
ave_pp = None
self._show_log_value_names()
if self._log_level > 2:
self._show_log_values_on_kstar(frequencies, gv, ave_pp, gp, q)
else:
self._show_log_values(frequencies, gv, ave_pp)
print("", end="", flush=True)
def _show_log_values(self, frequencies, gv, ave_pp):
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
for f, v, pp in zip(frequencies, gv, ave_pp):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e"
% (f, v[0], v[1], v[2], np.linalg.norm(v), pp)
)
else:
for f, v in zip(frequencies, gv):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f"
% (f, v[0], v[1], v[2], np.linalg.norm(v))
)
def _show_log_values_on_kstar(self, frequencies, gv, ave_pp, gp, q):
rotation_map = get_grid_points_by_rotations(gp, self._pp.bz_grid)
for i, j in enumerate(np.unique(rotation_map)):
for k, (rot, rot_c) in enumerate(
zip(self._point_operations, self._rotations_cartesian)
):
if rotation_map[k] != j:
continue
print(
" k*%-2d (%5.2f %5.2f %5.2f)" % ((i + 1,) + tuple(np.dot(rot, q)))
)
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
for f, v, pp in zip(frequencies, np.dot(rot_c, gv.T).T, ave_pp):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e"
% (f, v[0], v[1], v[2], np.linalg.norm(v), pp)
)
else:
for f, v in zip(frequencies, np.dot(rot_c, gv.T).T):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f"
% (f, v[0], v[1], v[2], np.linalg.norm(v))
)
print("")
def _show_log_value_names(self):
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
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:
pass
else:
text += " (dq=%3.1e)" % self._velocity_obj.q_length
print(text)
from phono3py.conductivity.base import ConductivityMixIn
from phono3py.conductivity.rta_base import ConductivityRTABase
from phono3py.phonon3.interaction import Interaction
class ConductivityRTA(ConductivityMixIn, ConductivityRTABase):
@ -649,473 +163,3 @@ class ConductivityRTA(ConductivityMixIn, ConductivityRTABase):
order="C",
dtype="double",
)
class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
"""Class of Wigner lattice thermal conductivity under RTA.
Authors
-------
Michele Simoncelli
"""
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,
):
"""Init method."""
self._cv = None
self._kappa_TOT_RTA = None
self._kappa_P_RTA = None
self._kappa_C = None
self._mode_kappa_P_RTA = None
self._mode_kappa_C = None
self._gv_operator = None
self._gv_operator_sum2 = None
super().__init__(
interaction,
grid_points=grid_points,
temperatures=temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
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,
pp_filename=pp_filename,
is_N_U=is_N_U,
is_gamma_detail=is_gamma_detail,
is_frequency_shift_by_bubble=is_frequency_shift_by_bubble,
log_level=log_level,
)
self._conversion_factor_WTE = get_conversion_factor_WTE(
self._pp.primitive.volume
)
def set_kappa_at_sigmas(self):
"""Calculate the Wigner thermal conductivity.
k_P + k_C using the scattering operator in the RTA approximation.
"""
num_band = len(self._pp.primitive) * 3
for i, _ in enumerate(self._grid_points):
cv = self._cv[:, i, :]
gp = self._grid_points[i]
frequencies = self._frequencies[gp]
# Kappa
for j in range(len(self._sigmas)):
for k in range(len(self._temperatures)):
g_sum = self._get_main_diagonal(
i, j, k
) # phonon HWHM: q-point, sigma, temperature
for s1 in range(num_band):
for s2 in range(num_band):
hbar_omega_eV_s1 = (
frequencies[s1] * THzToEv
) # hbar*omega=h*nu in eV
hbar_omega_eV_s2 = (
frequencies[s2] * THzToEv
) # hbar*omega=h*nu in eV
if (frequencies[s1] > self._pp.cutoff_frequency) and (
frequencies[s2] > self._pp.cutoff_frequency
):
hbar_gamma_eV_s1 = 2.0 * g_sum[s1] * THzToEv
hbar_gamma_eV_s2 = 2.0 * g_sum[s2] * THzToEv
#
lorentzian_divided_by_hbar = (
0.5 * (hbar_gamma_eV_s1 + hbar_gamma_eV_s2)
) / (
(hbar_omega_eV_s1 - hbar_omega_eV_s2) ** 2
+ 0.25
* ((hbar_gamma_eV_s1 + hbar_gamma_eV_s2) ** 2)
)
#
prefactor = (
0.25
* (hbar_omega_eV_s1 + hbar_omega_eV_s2)
* (
cv[k, s1] / hbar_omega_eV_s1
+ cv[k, 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])
* prefactor
* lorentzian_divided_by_hbar
* self._conversion_factor_WTE
).real
#
self._mode_kappa_P_RTA[j, k, i, s1] += (
0.5 * contribution
)
self._mode_kappa_P_RTA[j, k, i, s2] += (
0.5 * contribution
)
# prefactor 0.5 arises from the fact that degenerate
# modes have the same specific heat, hence they give
# the same contribution to the populations
# conductivity
else:
self._mode_kappa_C[j, k, i, s1, s2] += (
(self._gv_operator_sum2[i, s1, s2])
* prefactor
* lorentzian_divided_by_hbar
* self._conversion_factor_WTE
).real
elif s1 == s2:
self._num_ignored_phonon_modes[j, k] += 1
N = self._num_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 _allocate_values(self):
super()._allocate_values()
num_band0 = len(self._pp.band_indices)
num_grid_points = len(self._grid_points)
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(
(len(self._sigmas), num_temp, 6), order="C", dtype="double"
)
self._kappa_P_RTA = np.zeros(
(len(self._sigmas), num_temp, 6), order="C", dtype="double"
)
self._kappa_C = np.zeros(
(len(self._sigmas), num_temp, 6), order="C", dtype="double"
)
self._mode_kappa_P_RTA = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, nat3, 6),
order="C",
dtype="double",
)
# one more index because we have off-diagonal terms (second index not
# parallelized)
self._mode_kappa_C = np.zeros(
(
len(self._sigmas),
num_temp,
num_grid_points,
num_band0,
nat3,
6,
),
order="C",
dtype="double",
)
class ConductivityKuboRTA(ConductivityKuboMixIn, ConductivityRTABase):
"""Class of Kubo lattice thermal conductivity under 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,
):
"""Init method."""
self._cv_mat = None
self._gv_mat = None
self._kappa = None
super().__init__(
interaction,
grid_points=grid_points,
temperatures=temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
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,
pp_filename=pp_filename,
is_N_U=is_N_U,
is_gamma_detail=is_gamma_detail,
is_frequency_shift_by_bubble=is_frequency_shift_by_bubble,
log_level=log_level,
)
def set_kappa_at_sigmas(self):
"""Calculate kappa from ph-ph interaction results."""
for i_gp, _ in enumerate(self._grid_points):
frequencies = self._frequencies[self._grid_points[i_gp]]
for j in range(len(self._sigmas)):
for k in range(len(self._temperatures)):
g_sum = self._get_main_diagonal(i_gp, j, k)
for i_band, freq in enumerate(frequencies):
if freq < self._pp.cutoff_frequency:
self._num_ignored_phonon_modes[j, k] += 1
continue
self._set_kappa_at_sigmas(
j, k, i_gp, i_band, g_sum, frequencies
)
N = self._num_sampling_grid_points
self._kappa = self._mode_kappa_mat.sum(axis=2).sum(axis=2).sum(axis=2).real / N
def _set_kappa_at_sigmas(self, j, k, i_gp, i_band, g_sum, frequencies):
gvm_sum2 = self._gv_mat_sum2[i_gp]
cvm = self._cv_mat[k, i_gp]
for j_band, freq in enumerate(frequencies):
if freq < self._pp.cutoff_frequency:
return
g = g_sum[i_band] + g_sum[j_band]
for i_pair, _ in enumerate(
([0, 0], [1, 1], [2, 2], [1, 2], [0, 2], [0, 1])
):
old_settings = np.seterr(all="raise")
try:
self._mode_kappa_mat[j, k, i_gp, i_band, j_band, i_pair] = (
cvm[i_band, j_band]
* gvm_sum2[i_band, j_band, i_pair]
* g
/ ((frequencies[j_band] - frequencies[i_band]) ** 2 + g**2)
* self._conversion_factor
)
except FloatingPointError:
# supposed that g is almost 0 and |gv|=0
pass
except Exception:
gp = self._grid_points[i_gp]
print("=" * 26 + " Warning " + "=" * 26)
print(
" Unexpected physical condition of ph-ph "
"interaction calculation was found."
)
print(
" g=%f at gp=%d, band=%d, freq=%f, band=%d, freq=%f"
% (
g_sum[i_band],
gp,
i_band + 1,
frequencies[i_band],
j_band + 1,
frequencies[j_band],
)
)
print("=" * 61)
np.seterr(**old_settings)
def _allocate_values(self):
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_mat = np.zeros(
(num_temp, num_grid_points, num_band0, num_band), dtype="double", order="C"
)
self._gv_mat = np.zeros(
(num_grid_points, num_band0, num_band, 3),
dtype=self._complex_dtype,
order="C",
)
self._gv_mat_sum2 = np.zeros(
(num_grid_points, num_band0, num_band, 6),
dtype=self._complex_dtype,
order="C",
)
# kappa and mode_kappa_mat are accessed when all bands exist, i.e.,
# num_band0==num_band.
self._kappa = np.zeros(
(len(self._sigmas), num_temp, 6), dtype="double", order="C"
)
self._mode_kappa_mat = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band, num_band, 6),
dtype=self._complex_dtype,
order="C",
)
def get_thermal_conductivity_RTA(
interaction: Interaction,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
mass_variances=None,
grid_points=None,
is_isotope=False,
boundary_mfp=None, # in micrometer
use_ave_pp=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
is_N_U=False,
conductivity_type=None,
write_gamma=False,
read_gamma=False,
write_kappa=False,
write_pp=False,
read_pp=False,
write_gamma_detail=False,
compression="gzip",
input_filename=None,
output_filename=None,
log_level=0,
):
"""Run RTA thermal conductivity calculation."""
if temperatures is None:
_temperatures = np.arange(0, 1001, 10, dtype="double")
else:
_temperatures = temperatures
conductivity_RTA_class: Type[
Union[ConductivityRTA, ConductivityWignerRTA, ConductivityKuboRTA]
]
if conductivity_type == "wigner":
conductivity_RTA_class = ConductivityWignerRTA
elif conductivity_type == "kubo":
conductivity_RTA_class = ConductivityKuboRTA
else:
conductivity_RTA_class = ConductivityRTA
if log_level:
print(
"-------------------- Lattice thermal conductivity (RTA) "
"--------------------"
)
br = conductivity_RTA_class(
interaction,
grid_points=grid_points,
temperatures=_temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
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=write_pp,
pp_filename=input_filename,
is_N_U=is_N_U,
is_gamma_detail=write_gamma_detail,
log_level=log_level,
)
if read_gamma:
if not set_gamma_from_file(br, filename=input_filename):
print("Reading collisions failed.")
return False
for i in br:
if write_pp:
write_phph(
br, interaction, i, compression=compression, filename=output_filename
)
if write_gamma:
ConductivityRTAWriter.write_gamma(
br,
interaction,
i,
compression=compression,
filename=output_filename,
verbose=log_level,
)
if write_gamma_detail:
ConductivityRTAWriter.write_gamma_detail(
br,
interaction,
i,
compression=compression,
filename=output_filename,
verbose=log_level,
)
if grid_points is None and all_bands_exist(interaction):
br.set_kappa_at_sigmas()
if log_level:
if conductivity_type == "wigner":
ShowCalcProgress.kappa_Wigner_RTA(
cast(ConductivityWignerRTA, br), log_level
)
else:
ShowCalcProgress.kappa_RTA(cast(ConductivityRTA, br), log_level)
if write_kappa:
ConductivityRTAWriter.write_kappa(
br,
interaction.primitive.volume,
compression=compression,
filename=output_filename,
log_level=log_level,
)
return br

519
phono3py/conductivity/rta_base.py Normal file
View File

@ -0,0 +1,519 @@
"""Lattice thermal conductivity calculation base class with RTA."""
# Copyright (C) 2020 Atsushi Togo
# All rights reserved.
#
# This file is part of phono3py.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
#
# * Neither the name of the phonopy project nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import warnings
from abc import abstractmethod
import numpy as np
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.phonon3.imag_self_energy import ImagSelfEnergy, average_by_degeneracy
from phono3py.phonon3.interaction import Interaction
class ConductivityRTABase(ConductivityBase):
"""Base class of ConductivityRTA*.
This is a base class of RTA classes.
"""
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,
):
"""Init method."""
self._is_N_U = is_N_U
self._is_gamma_detail = is_gamma_detail
self._is_frequency_shift_by_bubble = is_frequency_shift_by_bubble
self._gamma_N = None
self._gamma_U = None
self._gamma_detail_at_q = None
self._use_ave_pp = use_ave_pp
self._use_const_ave_pp = None
self._num_ignored_phonon_modes = None
super().__init__(
interaction,
grid_points=grid_points,
temperatures=temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
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,
)
self._use_const_ave_pp = self._pp.constant_averaged_interaction
self._read_pp = read_pp
self._store_pp = store_pp
self._pp_filename = pp_filename
if self._temperatures is not None:
self._allocate_values()
self._collision = ImagSelfEnergy(
self._pp, with_detail=(self._is_gamma_detail or self._is_N_U)
)
def get_gamma_N_U(self):
"""Return N and U parts of gamma."""
return (self._gamma_N, self._gamma_U)
def set_gamma_N_U(self, gamma_N, gamma_U):
"""Set N and U parts of gamma."""
self._gamma_N = gamma_N
self._gamma_U = gamma_U
def get_gamma_detail_at_q(self):
"""Return contribution of each triplet to gamma at current q-point."""
return self._gamma_detail_at_q
def get_number_of_ignored_phonon_modes(self):
"""Return number of ignored phonon modes."""
warnings.warn(
"Use attribute, number_of_ignored_phonon_modes",
DeprecationWarning,
stacklevel=2,
)
return self.number_of_ignored_phonon_modes
@property
def number_of_ignored_phonon_modes(self):
"""Return number of ignored phonon modes."""
return self._num_ignored_phonon_modes
def set_averaged_pp_interaction(self, ave_pp):
"""Set averaged ph-ph interaction."""
self._averaged_pp_interaction = ave_pp
@abstractmethod
def set_kappa_at_sigmas(self):
"""Must be implemented in the inherited class."""
raise NotImplementedError()
def _allocate_values(self):
num_band0 = len(self._pp.band_indices)
num_grid_points = len(self._grid_points)
num_temp = len(self._temperatures)
if not self._read_gamma:
self._gamma = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band0),
order="C",
dtype="double",
)
if self._is_gamma_detail or self._is_N_U:
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),
order="C",
dtype="double",
)
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
self._averaged_pp_interaction = np.zeros(
(num_grid_points, num_band0), order="C", dtype="double"
)
self._num_ignored_phonon_modes = np.zeros(
(len(self._sigmas), num_temp), order="C", dtype="intc"
)
def _run_at_grid_point(self):
i_gp = self._grid_point_count
self._show_log_header(i_gp)
grid_point = self._grid_points[i_gp]
self._set_cv(i_gp, i_gp)
self._set_velocities(i_gp, i_gp)
if self._read_gamma:
if self._use_ave_pp:
self._collision.set_grid_point(grid_point)
self._set_gamma_at_sigmas(i_gp)
else:
self._collision.set_grid_point(grid_point)
num_triplets = len(self._pp.get_triplets_at_q()[0])
if self._log_level:
print("Number of triplets: %d" % num_triplets)
if (
self._is_full_pp
or self._read_pp
or self._store_pp
or self._use_ave_pp
or self._use_const_ave_pp
or self._is_gamma_detail
):
self._set_gamma_at_sigmas(i_gp)
else: # can save memory space
self._set_gamma_at_sigmas_lowmem(i_gp)
if self._is_isotope and not self._read_gamma_iso:
gamma_iso = self._get_gamma_isotope_at_sigmas(i_gp)
self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._pp.band_indices]
if self._log_level:
self._show_log(i_gp)
def _set_gamma_at_sigmas(self, i):
for j, sigma in enumerate(self._sigmas):
self._collision.set_sigma(sigma, sigma_cutoff=self._sigma_cutoff)
self._collision.run_integration_weights()
if self._log_level:
text = "Collisions will be calculated with "
if sigma is None:
text += "tetrahedron method."
else:
text += "sigma=%s" % sigma
if self._sigma_cutoff is None:
text += "."
else:
text += "(%4.2f SD)." % self._sigma_cutoff
print(text)
if self._read_pp:
pp, _g_zero = read_pp_from_hdf5(
self._pp.mesh_numbers,
grid_point=self._grid_points[i],
sigma=sigma,
sigma_cutoff=self._sigma_cutoff,
filename=self._pp_filename,
verbose=(self._log_level > 0),
)
_, g_zero = self._collision.get_integration_weights()
if self._log_level:
if len(self._sigmas) > 1:
print(
"Multiple sigmas or mixing smearing and "
"tetrahedron method is not supported."
)
if _g_zero is not None and (_g_zero != g_zero).any():
raise ValueError("Inconsistency found in g_zero.")
self._collision.set_interaction_strength(pp)
elif self._use_ave_pp:
self._collision.set_averaged_pp_interaction(
self._averaged_pp_interaction[i]
)
elif self._use_const_ave_pp:
if self._log_level:
print(
"Constant ph-ph interaction of %6.3e is used."
% self._pp.constant_averaged_interaction
)
self._collision.run_interaction()
self._averaged_pp_interaction[i] = self._pp.averaged_interaction
elif j != 0 and (self._is_full_pp or self._sigma_cutoff is None):
if self._log_level:
print("Existing ph-ph interaction is used.")
else:
if self._log_level:
print("Calculating ph-ph interaction...")
self._collision.run_interaction(is_full_pp=self._is_full_pp)
if self._is_full_pp:
self._averaged_pp_interaction[i] = self._pp.averaged_interaction
# Number of triplets depends on q-point.
# So this is allocated each time.
if self._is_gamma_detail:
num_temp = len(self._temperatures)
self._gamma_detail_at_q = np.empty(
((num_temp,) + self._pp.get_interaction_strength().shape),
dtype="double",
order="C",
)
self._gamma_detail_at_q[:] = 0
if self._log_level:
print("Calculating collisions at temperatures...")
for k, t in enumerate(self._temperatures):
self._collision.temperature = t
self._collision.run()
self._gamma[j, k, i] = self._collision.imag_self_energy
if self._is_N_U or self._is_gamma_detail:
g_N, g_U = self._collision.get_imag_self_energy_N_and_U()
self._gamma_N[j, k, i] = g_N
self._gamma_U[j, k, i] = g_U
if self._is_gamma_detail:
self._gamma_detail_at_q[k] = (
self._collision.get_detailed_imag_self_energy()
)
def _set_gamma_at_sigmas_lowmem(self, i):
"""Calculate gamma without storing ph-ph interaction strength.
`svecs` and `multi` below must not be simply replaced by
`self._pp.primitive.get_smallest_vectors()` because they must be in
dense format as always so in Interaction class instance.
`p2s`, `s2p`, and `masses` have to be also given from Interaction
class instance.
"""
num_band = len(self._pp.primitive) * 3
band_indices = self._pp.band_indices
(
svecs,
multi,
p2s,
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:
tetrahedra = get_tetrahedra_relative_grid_address(
self._pp.bz_grid.microzone_lattice
)
# It is assumed that self._sigmas = [None].
temperatures_THz = np.array(
self._temperatures * get_physical_units().KB / get_physical_units().THzToEv,
dtype="double",
)
for j, sigma in enumerate(self._sigmas):
self._collision.set_sigma(sigma)
if self._is_N_U:
collisions = np.zeros(
(2, len(self._temperatures), len(band_indices)),
dtype="double",
order="C",
)
else:
collisions = np.zeros(
(len(self._temperatures), len(band_indices)),
dtype="double",
order="C",
)
import phono3py._phono3py as phono3c
# True: OpenMP over triplets
# False: OpenMP over bands
if self._pp.openmp_per_triplets is None:
if len(triplets_at_q) > num_band:
openmp_per_triplets = True
else:
openmp_per_triplets = False
else:
openmp_per_triplets = self._pp.openmp_per_triplets
if sigma is None:
phono3c.pp_collision(
collisions,
np.array(
np.dot(tetrahedra, self._pp.bz_grid.P.T),
dtype="int64",
order="C",
),
self._frequencies,
self._eigenvectors,
triplets_at_q,
weights_at_q,
self._pp.bz_grid.addresses,
self._pp.bz_grid.gp_map,
self._pp.bz_grid.store_dense_gp_map * 1 + 1,
self._pp.bz_grid.D_diag,
self._pp.bz_grid.Q,
fc3,
svecs,
multi,
masses,
p2s,
s2p,
band_indices,
temperatures_THz,
self._is_N_U * 1,
self._pp.symmetrize_fc3q * 1,
self._pp.make_r0_average * 1,
self._pp.all_shortest,
self._pp.cutoff_frequency,
openmp_per_triplets * 1,
)
else:
if self._sigma_cutoff is None:
sigma_cutoff = -1
else:
sigma_cutoff = float(self._sigma_cutoff)
phono3c.pp_collision_with_sigma(
collisions,
sigma,
sigma_cutoff,
self._frequencies,
self._eigenvectors,
triplets_at_q,
weights_at_q,
self._pp.bz_grid.addresses,
self._pp.bz_grid.D_diag,
self._pp.bz_grid.Q,
fc3,
svecs,
multi,
masses,
p2s,
s2p,
band_indices,
temperatures_THz,
self._is_N_U * 1,
self._pp.symmetrize_fc3q * 1,
self._pp.make_r0_average * 1,
self._pp.all_shortest,
self._pp.cutoff_frequency,
openmp_per_triplets * 1,
)
col_unit_conv = self._collision.unit_conversion_factor
pp_unit_conv = self._pp.unit_conversion_factor
if self._is_N_U:
col = collisions.sum(axis=0)
col_N = collisions[0]
col_U = collisions[1]
else:
col = collisions
for k in range(len(self._temperatures)):
self._gamma[j, k, i, :] = average_by_degeneracy(
col[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]],
)
if self._is_N_U:
self._gamma_N[j, k, i, :] = average_by_degeneracy(
col_N[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]],
)
self._gamma_U[j, k, i, :] = average_by_degeneracy(
col_U[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]],
)
def _show_log(self, i_gp):
q = self._get_qpoint_from_gp_index(i_gp)
gp = self._grid_points[i_gp]
frequencies = self._frequencies[gp][self._pp.band_indices]
gv = self._gv[i_gp]
if self._averaged_pp_interaction is not None:
ave_pp = self._averaged_pp_interaction[i_gp]
else:
ave_pp = None
self._show_log_value_names()
if self._log_level > 2:
self._show_log_values_on_kstar(frequencies, gv, ave_pp, gp, q)
else:
self._show_log_values(frequencies, gv, ave_pp)
print("", end="", flush=True)
def _show_log_values(self, frequencies, gv, ave_pp):
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
for f, v, pp in zip(frequencies, gv, ave_pp):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e"
% (f, v[0], v[1], v[2], np.linalg.norm(v), pp)
)
else:
for f, v in zip(frequencies, gv):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f"
% (f, v[0], v[1], v[2], np.linalg.norm(v))
)
def _show_log_values_on_kstar(self, frequencies, gv, ave_pp, gp, q):
rotation_map = get_grid_points_by_rotations(gp, self._pp.bz_grid)
for i, j in enumerate(np.unique(rotation_map)):
for k, (rot, rot_c) in enumerate(
zip(self._point_operations, self._rotations_cartesian)
):
if rotation_map[k] != j:
continue
print(
" k*%-2d (%5.2f %5.2f %5.2f)" % ((i + 1,) + tuple(np.dot(rot, q)))
)
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
for f, v, pp in zip(frequencies, np.dot(rot_c, gv.T).T, ave_pp):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e"
% (f, v[0], v[1], v[2], np.linalg.norm(v), pp)
)
else:
for f, v in zip(frequencies, np.dot(rot_c, gv.T).T):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f"
% (f, v[0], v[1], v[2], np.linalg.norm(v))
)
print("")
def _show_log_value_names(self):
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
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:
pass
else:
text += " (dq=%3.1e)" % self._velocity_obj.q_length
print(text)

1144
phono3py/conductivity/utils.py
View File

File diff suppressed because it is too large Load Diff

13
phono3py/conductivity/wigner.py → phono3py/conductivity/wigner_base.py
View File

@ -38,7 +38,7 @@ import textwrap
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.units import EV, Angstrom, Hbar, THz
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.base import HeatCapacityMixIn
from phono3py.phonon.grid import get_grid_points_by_rotations
@ -215,8 +215,9 @@ class ConductivityWignerMixIn(HeatCapacityMixIn):
def get_conversion_factor_WTE(volume):
"""Return conversion factor of thermal conductivity."""
return (
(THz * Angstrom) ** 2 # ----> group velocity
* EV # ----> specific heat is in eV/
* Hbar # ----> transform lorentzian_div_hbar from eV^-1 to s
/ (volume * Angstrom**3)
) # ----> unit cell volume
(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

328
phono3py/conductivity/wigner_direct_solution.py Normal file
View File

@ -0,0 +1,328 @@
"""Wigner thermal conductivity direct solution class."""
# Copyright (C) 2022 Michele Simoncelli
# All rights reserved.
#
# This file is part of phono3py.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
#
# * Neither the name of the phonopy project nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.physical_units import get_physical_units
from phono3py.conductivity.direct_solution import (
ConductivityLBTEBase,
diagonalize_collision_matrix,
)
from phono3py.conductivity.wigner_base import (
ConductivityWignerMixIn,
get_conversion_factor_WTE,
)
from phono3py.phonon3.interaction import Interaction
class ConductivityWignerLBTE(ConductivityWignerMixIn, ConductivityLBTEBase):
"""Class of Wigner lattice thermal conductivity under direct-solution.
Authors
-------
Michele Simoncelli
"""
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",
):
"""Init method."""
self._kappa_TOT_exact = None
self._kappa_TOT_RTA = None
self._kappa_P_exact = None
self._mode_kappa_P_exact = None
self._kappa_P_RTA = None
self._mode_kappa_P_RTA = None
self._kappa_C = None
self._mode_kappa_C = None
self._gv_operator = None
self._gv_operator_sum2 = None
super().__init__(
interaction,
grid_points=grid_points,
temperatures=temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
boundary_mfp=boundary_mfp,
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,
pinv_cutoff=pinv_cutoff,
pinv_solver=pinv_solver,
pinv_method=pinv_method,
log_level=log_level,
lang=lang,
)
self._conversion_factor_WTE = get_conversion_factor_WTE(
self._pp.primitive.volume
)
@property
def kappa_TOT_exact(self):
"""Return kappa."""
return self._kappa_TOT_exact
@property
def kappa_P_exact(self):
"""Return kappa."""
return self._kappa_P_exact
@property
def mode_kappa_P_exact(self):
"""Return mode_kappa."""
return self._mode_kappa_P_exact
def _allocate_local_values(self, num_grid_points):
"""Allocate grid point local arrays."""
num_band0 = len(self._pp.band_indices)
num_temp = len(self._temperatures)
super()._allocate_local_values(num_grid_points)
nat3 = len(self._pp.primitive) * 3
self._kappa_TOT_RTA = np.zeros(
(len(self._sigmas), num_temp, 6), dtype="double", order="C"
)
self._kappa_TOT_exact = np.zeros(
(len(self._sigmas), num_temp, 6), dtype="double", order="C"
)
self._kappa_P_exact = np.zeros(
(len(self._sigmas), num_temp, 6), dtype="double", order="C"
)
self._kappa_P_RTA = np.zeros(
(len(self._sigmas), num_temp, 6), dtype="double", order="C"
)
self._kappa_C = np.zeros(
(len(self._sigmas), num_temp, 6), dtype="double", order="C"
)
self._mode_kappa_P_exact = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band0, 6), dtype="double"
)
self._mode_kappa_P_RTA = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band0, 6), dtype="double"
)
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(
(
len(self._sigmas),
num_temp,
num_grid_points,
num_band0,
nat3,
6,
),
order="C",
dtype=complex_dtype,
)
def _set_kappa_at_sigmas(self, weights):
"""Calculate thermal conductivity from collision matrix."""
for j, sigma in enumerate(self._sigmas):
if self._log_level:
text = "----------- Thermal conductivity (W/m-k) "
if sigma:
text += "for sigma=%s -----------" % sigma
else:
text += "with tetrahedron method -----------"
print(text, flush=True)
for k, t in enumerate(self._temperatures):
if t > 0:
self._set_kappa_RTA(j, k, weights)
w = diagonalize_collision_matrix(
self._collision_matrix,
i_sigma=j,
i_temp=k,
pinv_solver=self._pinv_solver,
log_level=self._log_level,
)
self._collision_eigenvalues[j, k] = w
self._set_kappa(j, k, weights)
if self._log_level:
print(
("#%6s " + " %-10s" * 6)
% (
" \t\t T(K)",
"xx",
"yy",
"zz",
"yz",
"xz",
"xy",
)
)
print(
"K_P_exact\t\t"
+ ("%7.1f " + " %10.3f" * 6)
% ((t,) + tuple(self._kappa_P_exact[j, k]))
)
print(
"(K_P_RTA)\t\t"
+ ("%7.1f " + " %10.3f" * 6)
% ((t,) + tuple(self._kappa_P_RTA[j, k]))
)
print(
"K_C \t\t"
+ ("%7.1f " + " %10.3f" * 6)
% ((t,) + tuple(self._kappa_C[j, k].real))
)
print(" ")
print(
"K_TOT=K_P_exact+K_C\t"
+ ("%7.1f " + " %10.3f" * 6)
% (
(t,)
+ tuple(self._kappa_C[j, k] + self._kappa_P_exact[j, k])
)
)
print("-" * 76, flush=True)
if self._log_level:
print("", flush=True)
def _set_kappa(self, i_sigma, i_temp, weights):
if self._is_reducible_collision_matrix:
self._set_kappa_reducible_colmat(
self._kappa_P_exact, self._mode_kappa_P_exact, i_sigma, i_temp, weights
)
else:
self._set_kappa_ir_colmat(
self._kappa_P_exact, self._mode_kappa_P_exact, i_sigma, i_temp, weights
)
def _set_kappa_RTA(self, i_sigma, i_temp, weights):
if self._is_reducible_collision_matrix:
self._set_kappa_RTA_reducible_colmat(
self._kappa_P_RTA, self._mode_kappa_P_RTA, i_sigma, i_temp, weights
)
print(
" WARNING: Coherences conductivity not implemented for "
"is_reducible_collision_matrix=True "
)
else:
self._set_kappa_RTA_ir_colmat(
self._kappa_P_RTA, self._mode_kappa_P_RTA, i_sigma, i_temp, weights
)
self._set_kappa_C_ir_colmat(i_sigma, i_temp)
def _set_kappa_C_ir_colmat(self, i_sigma, i_temp):
"""Calculate coherence term of the thermal conductivity."""
N = self._num_sampling_grid_points
num_band = len(self._pp.primitive) * 3
# num_ir_grid_points = len(self._ir_grid_points)
THzToEv = get_physical_units().THzToEv
for i, gp in enumerate(self._ir_grid_points):
# 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, :]
for s1 in range(num_band):
for s2 in range(num_band):
hbar_omega_eV_s1 = (
frequencies[s1] * THzToEv
) # hbar*omega=h*nu in eV
hbar_omega_eV_s2 = (
frequencies[s2] * THzToEv
) # hbar*omega=h*nu in eV
if (
(frequencies[s1] > self._pp.cutoff_frequency)
and (frequencies[s2] > self._pp.cutoff_frequency)
) and np.abs(frequencies[s1] - frequencies[s2]) > 1e-4:
# frequencies must be non-degenerate to contribute to k_C,
# otherwise they contribute to k_P
hbar_gamma_eV_s1 = g[s1] * THzToEv
hbar_gamma_eV_s2 = g[s2] * THzToEv
lorentzian_divided_by_hbar = (
0.5 * (hbar_gamma_eV_s1 + hbar_gamma_eV_s2)
) / (
(hbar_omega_eV_s1 - hbar_omega_eV_s2) ** 2
+ 0.25 * ((hbar_gamma_eV_s1 + hbar_gamma_eV_s2) ** 2)
)
#
prefactor = (
0.25
* (hbar_omega_eV_s1 + hbar_omega_eV_s2)
* (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])
* prefactor
* lorentzian_divided_by_hbar
* self._conversion_factor_WTE
)
self._kappa_C[i_sigma, i_temp] = (
self._mode_kappa_C[i_sigma, i_temp].sum(axis=0).sum(axis=0).sum(axis=0) / N
).real

250
phono3py/conductivity/wigner_rta.py Normal file
View File

@ -0,0 +1,250 @@
"""Wigner thermal conductivity RTA class."""
# Copyright (C) 2022 Michele Simoncelli
# All rights reserved.
#
# This file is part of phono3py.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in
# the documentation and/or other materials provided with the
# distribution.
#
# * Neither the name of the phonopy project nor the names of its
# contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
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,
get_conversion_factor_WTE,
)
from phono3py.phonon3.interaction import Interaction
class ConductivityWignerRTA(ConductivityWignerMixIn, ConductivityRTABase):
"""Class of Wigner lattice thermal conductivity under RTA.
Authors
-------
Michele Simoncelli
"""
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,
):
"""Init method."""
self._cv = None
self._kappa_TOT_RTA = None
self._kappa_P_RTA = None
self._kappa_C = None
self._mode_kappa_P_RTA = None
self._mode_kappa_C = None
self._gv_operator = None
self._gv_operator_sum2 = None
super().__init__(
interaction,
grid_points=grid_points,
temperatures=temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
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,
pp_filename=pp_filename,
is_N_U=is_N_U,
is_gamma_detail=is_gamma_detail,
is_frequency_shift_by_bubble=is_frequency_shift_by_bubble,
log_level=log_level,
)
self._conversion_factor_WTE = get_conversion_factor_WTE(
self._pp.primitive.volume
)
def set_kappa_at_sigmas(self):
"""Calculate the Wigner thermal conductivity.
k_P + k_C using the scattering operator in the RTA approximation.
"""
num_band = len(self._pp.primitive) * 3
THzToEv = get_physical_units().THzToEv
for i, _ in enumerate(self._grid_points):
cv = self._cv[:, i, :]
gp = self._grid_points[i]
frequencies = self._frequencies[gp]
# Kappa
for j in range(len(self._sigmas)):
for k in range(len(self._temperatures)):
g_sum = self._get_main_diagonal(
i, j, k
) # phonon HWHM: q-point, sigma, temperature
for s1 in range(num_band):
for s2 in range(num_band):
hbar_omega_eV_s1 = (
frequencies[s1] * THzToEv
) # hbar*omega=h*nu in eV
hbar_omega_eV_s2 = (
frequencies[s2] * THzToEv
) # hbar*omega=h*nu in eV
if (frequencies[s1] > self._pp.cutoff_frequency) and (
frequencies[s2] > self._pp.cutoff_frequency
):
hbar_gamma_eV_s1 = 2.0 * g_sum[s1] * THzToEv
hbar_gamma_eV_s2 = 2.0 * g_sum[s2] * THzToEv
#
lorentzian_divided_by_hbar = (
0.5 * (hbar_gamma_eV_s1 + hbar_gamma_eV_s2)
) / (
(hbar_omega_eV_s1 - hbar_omega_eV_s2) ** 2
+ 0.25
* ((hbar_gamma_eV_s1 + hbar_gamma_eV_s2) ** 2)
)
#
prefactor = (
0.25
* (hbar_omega_eV_s1 + hbar_omega_eV_s2)
* (
cv[k, s1] / hbar_omega_eV_s1
+ cv[k, 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])
* prefactor
* lorentzian_divided_by_hbar
* self._conversion_factor_WTE
).real
#
self._mode_kappa_P_RTA[j, k, i, s1] += (
0.5 * contribution
)
self._mode_kappa_P_RTA[j, k, i, s2] += (
0.5 * contribution
)
# prefactor 0.5 arises from the fact that degenerate
# modes have the same specific heat, hence they give
# the same contribution to the populations
# conductivity
else:
self._mode_kappa_C[j, k, i, s1, s2] += (
(self._gv_operator_sum2[i, s1, s2])
* prefactor
* lorentzian_divided_by_hbar
* self._conversion_factor_WTE
).real
elif s1 == s2:
self._num_ignored_phonon_modes[j, k] += 1
N = self._num_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 _allocate_values(self):
super()._allocate_values()
num_band0 = len(self._pp.band_indices)
num_grid_points = len(self._grid_points)
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(
(len(self._sigmas), num_temp, 6), order="C", dtype="double"
)
self._kappa_P_RTA = np.zeros(
(len(self._sigmas), num_temp, 6), order="C", dtype="double"
)
self._kappa_C = np.zeros(
(len(self._sigmas), num_temp, 6), order="C", dtype="double"
)
self._mode_kappa_P_RTA = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, nat3, 6),
order="C",
dtype="double",
)
# one more index because we have off-diagonal terms (second index not
# parallelized)
self._mode_kappa_C = np.zeros(
(
len(self._sigmas),
num_temp,
num_grid_points,
num_band0,
nat3,
6,
),
order="C",
dtype="double",
)

4
phono3py/cui/create_force_constants.py
View File

@ -50,7 +50,7 @@ from phonopy.harmonic.force_constants import (
symmetrize_compact_force_constants,
symmetrize_force_constants,
)
from phonopy.interface.calculator import get_default_physical_units
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
@ -254,7 +254,7 @@ def parse_forces(
if dataset:
filename_read_from = phono3py_yaml_filename
physical_units = get_default_physical_units(calculator)
physical_units = get_calculator_physical_units(calculator)
# Forces are not yet found in dataset. Then try to read from FORCES_FC3 or
# FORCES_FC2.

14
phono3py/cui/create_force_sets.py
View File

@ -477,9 +477,10 @@ def _set_forces_and_nac_params(
ph3py_yaml.dataset["forces"] = np.array(
calc_dataset_fc3["forces"], dtype="double", order="C"
)
ph3py_yaml.dataset["supercell_energies"] = np.array(
calc_dataset_fc3["supercell_energies"], dtype="double"
)
if "supercell_energies" in calc_dataset_fc3:
ph3py_yaml.dataset["supercell_energies"] = np.array(
calc_dataset_fc3["supercell_energies"], dtype="double"
)
if len(ph3py_yaml.dataset["forces"]) != len(
ph3py_yaml.dataset["displacements"]
):
@ -499,9 +500,10 @@ def _set_forces_and_nac_params(
ph3py_yaml.phonon_dataset["forces"] = np.array(
calc_dataset_fc2["forces"], dtype="double", order="C"
)
ph3py_yaml.phonon_dataset["supercell_energies"] = np.array(
calc_dataset_fc2["supercell_energies"], dtype="double"
)
if "supercell_energies" in calc_dataset_fc2:
ph3py_yaml.phonon_dataset["supercell_energies"] = np.array(
calc_dataset_fc2["supercell_energies"], dtype="double"
)
if len(ph3py_yaml.phonon_dataset["forces"]) != len(
ph3py_yaml.phonon_dataset["displacements"]
):

8
phono3py/cui/load.py
View File

@ -43,10 +43,10 @@ from typing import Optional, Union
import numpy as np
import phonopy.cui.load_helper as load_helper
from phonopy.harmonic.force_constants import show_drift_force_constants
from phonopy.interface.calculator import get_default_physical_units
from phonopy.interface.calculator import get_calculator_physical_units
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import determinant
from phonopy.units import VaspToTHz
from phono3py import Phono3py
from phono3py.cui.create_force_constants import (
@ -310,12 +310,12 @@ def load(
_nac_params = None
# Convert distance unit of unit cell to Angstrom
physical_units = get_default_physical_units(_calculator)
physical_units = get_calculator_physical_units(_calculator)
factor_to_A = physical_units["distance_to_A"]
cell.cell = cell.cell * factor_to_A
if factor is None:
_factor = VaspToTHz
_factor = get_physical_units().DefaultToTHz
else:
_factor = factor
ph3py = Phono3py(

16
phono3py/cui/phono3py_script.py
View File

@ -58,10 +58,10 @@ from phonopy.cui.phonopy_script import (
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_default_physical_units
from phonopy.interface.calculator import get_calculator_physical_units
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
from phonopy.units import Bohr, Hartree, VaspToTHz
from phono3py import Phono3py, Phono3pyIsotope, Phono3pyJointDos
from phono3py.cui.create_force_constants import (
@ -159,7 +159,7 @@ def finalize_phono3py(
else:
yaml_filename = filename
_physical_units = get_default_physical_units(phono3py.calculator)
_physical_units = get_calculator_physical_units(phono3py.calculator)
ph3py_yaml = Phono3pyYaml(
configuration=confs_dict,
@ -390,7 +390,7 @@ def get_default_values(settings):
temperatures = settings.temperatures # For others
if settings.frequency_conversion_factor is None:
frequency_factor_to_THz = VaspToTHz
frequency_factor_to_THz = get_physical_units().DefaultToTHz
else:
frequency_factor_to_THz = settings.frequency_conversion_factor
@ -475,7 +475,7 @@ def init_phono3py(
settings, cell_info, interface_mode, symprec, log_level
) -> tuple[Phono3py, dict]:
"""Initialize phono3py and update settings by default values."""
physical_units = get_default_physical_units(interface_mode)
physical_units = get_calculator_physical_units(interface_mode)
distance_to_A = physical_units["distance_to_A"]
# Change unit of lattice parameters to angstrom
@ -566,7 +566,7 @@ def create_supercells_with_displacements(
cell_info["phonopy_yaml"],
unitcell_filename,
log_level,
nac_factor=Hartree * Bohr,
nac_factor=get_physical_units().Hartree * get_physical_units().Bohr,
load_phonopy_yaml=load_phono3py_yaml,
)
@ -730,7 +730,7 @@ def run_gruneisen_then_exit(phono3py, settings, output_filename, log_level):
nac_params=phono3py.nac_params,
nac_q_direction=settings.nac_q_direction,
ion_clamped=settings.ion_clamped,
factor=VaspToTHz,
factor=get_physical_units().DefaultToTHz,
symprec=phono3py.symmetry.tolerance,
output_filename=output_filename,
log_level=log_level,
@ -1149,7 +1149,7 @@ def main(**argparse_control):
cell_info["phonopy_yaml"],
unitcell_filename,
log_level,
nac_factor=Hartree * Bohr,
nac_factor=get_physical_units().Hartree * get_physical_units().Bohr,
load_phonopy_yaml=load_phono3py_yaml,
)

9
phono3py/other/isotope.py
View File

@ -41,12 +41,12 @@ from typing import Optional, Union
import numpy as np
from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
from phonopy.phonon.tetrahedron_mesh import get_tetrahedra_frequencies
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.atoms import isotope_data as phonopy_isotope_data
from phonopy.structure.cells import Primitive
from phonopy.structure.symmetry import Symmetry
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phonopy.units import VaspToTHz
from phono3py.other.tetrahedron_method import (
get_integration_weights,
@ -100,7 +100,7 @@ class Isotope:
band_indices=None,
sigma=None,
bz_grid=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
use_grg=False,
symprec=1e-5,
cutoff_frequency=None,
@ -122,7 +122,10 @@ class Isotope:
self._cutoff_frequency = 0
else:
self._cutoff_frequency = cutoff_frequency
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().DefaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._lapack_zheev_uplo = lapack_zheev_uplo
self._nac_q_direction = None

14
phono3py/phonon/func.py
View File

@ -35,7 +35,7 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.units import AMU, EV, Angstrom, Hbar, Kb, THz, THzToEv
from phonopy.physical_units import get_physical_units
def gaussian(x, sigma):
@ -62,7 +62,9 @@ def bose_einstein(x, T):
Temperature in K
"""
return 1.0 / (np.exp(THzToEv * x / (Kb * T)) - 1)
return 1.0 / (
np.exp(get_physical_units().THzToEv * x / (get_physical_units().KB * T)) - 1
)
def sigma_squared(x, T):
@ -96,7 +98,13 @@ def sigma_squared(x, T):
n = bose_einstein(x, T)
# factor=1.0107576777968994
factor = Hbar * EV / (2 * np.pi * THz) / AMU / Angstrom**2
factor = (
get_physical_units().Hbar
* get_physical_units().EV
/ (2 * np.pi * get_physical_units().THz)
/ get_physical_units().AMU
/ get_physical_units().Angstrom ** 2
)
#########################
np.seterr(**old_settings)

3
phono3py/phonon/group_velocity_matrix.py
View File

@ -37,7 +37,6 @@
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.units import VaspToTHz
from phonopy.utils import similarity_transformation
@ -57,7 +56,7 @@ class GroupVelocityMatrix(GroupVelocity):
dynamical_matrix,
q_length=None,
symmetry=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
cutoff_frequency=1e-4,
):
"""Init method.

7
phono3py/phonon/heat_capacity_matrix.py
View File

@ -36,7 +36,7 @@
import numpy as np
from phonopy.phonon.thermal_properties import mode_cv
from phonopy.units import Kb
from phonopy.physical_units import get_physical_units
def mode_cv_matrix(temp, freqs, cutoff=1e-4):
@ -67,7 +67,8 @@ def mode_cv_matrix(temp, freqs, cutoff=1e-4):
shape=(num_band, num_band), dtype='double', order='C'.
"""
x = freqs / Kb / temp
KB = get_physical_units().KB
x = freqs / KB / temp
shape = (len(freqs), len(freqs))
cvm = np.zeros(shape, dtype="double", order="C")
for i, j in np.ndindex(shape):
@ -77,5 +78,5 @@ def mode_cv_matrix(temp, freqs, cutoff=1e-4):
sub = x[i] - x[j]
add = x[i] + x[j]
n_inv = np.exp([x[i], x[j], sub]) - 1
cvm[i, j] = Kb * n_inv[2] / sub * (add / 2) ** 2 / n_inv[0] * (1 / n_inv[1] + 1)
cvm[i, j] = KB * n_inv[2] / sub * (add / 2) ** 2 / n_inv[0] * (1 / n_inv[1] + 1)
return cvm

13
phono3py/phonon/solver.py
View File

@ -35,8 +35,8 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import sparse_to_dense_svecs
from phonopy.units import VaspToTHz
def run_phonon_solver_c(
@ -47,7 +47,7 @@ def run_phonon_solver_c(
grid_points,
grid_address,
QDinv,
frequency_conversion_factor=VaspToTHz,
frequency_conversion_factor=None,
nac_q_direction=None, # in reduced coordinates
lapack_zheev_uplo="L",
):
@ -84,6 +84,11 @@ def run_phonon_solver_c(
import phono3py._phono3py as phono3c
import phono3py._phononcalc as phononcalc
if frequency_conversion_factor is None:
_frequency_conversion_factor = get_physical_units().DefaultToTHz
else:
_frequency_conversion_factor = frequency_conversion_factor
(
svecs,
multi,
@ -154,7 +159,7 @@ def run_phonon_solver_c(
masses,
fc_p2s,
fc_s2p,
frequency_conversion_factor,
_frequency_conversion_factor,
born,
dielectric,
rec_lattice,
@ -179,7 +184,7 @@ def run_phonon_solver_c(
frequencies[phonon_undone] = (
np.sign(frequencies[phonon_undone])
* np.sqrt(np.abs(frequencies[phonon_undone]))
* frequency_conversion_factor
* _frequency_conversion_factor
)

9
phono3py/phonon/velocity_operator.py
View File

@ -36,7 +36,7 @@
import numpy as np
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.units import VaspToTHz
from phonopy.physical_units import get_physical_units
class VelocityOperator(GroupVelocity):
@ -47,7 +47,7 @@ class VelocityOperator(GroupVelocity):
dynamical_matrix,
q_length=None,
symmetry=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
cutoff_frequency=1e-4,
):
"""Init method.
@ -75,7 +75,10 @@ class VelocityOperator(GroupVelocity):
if self._q_length is None:
self._q_length = 5e-6
self._symmetry = symmetry
self._factor = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._factor = get_physical_units().DefaultToTHz
else:
self._factor = frequency_factor_to_THz
self._cutoff_frequency = cutoff_frequency
self._directions = np.array(

12
phono3py/phonon3/collision_matrix.py
View File

@ -37,7 +37,7 @@
from typing import Optional
import numpy as np
from phonopy.units import Kb, THzToEv
from phonopy.physical_units import get_physical_units
from phono3py.phonon3.imag_self_energy import ImagSelfEnergy
from phono3py.phonon3.interaction import Interaction
@ -164,7 +164,7 @@ class CollisionMatrix(ImagSelfEnergy):
self._ir_map_at_q,
self._rot_grid_points, # in GRGrid
self._rotations_cartesian,
self._temperature,
self._temperature * get_physical_units().KB / get_physical_units().THzToEv,
self._unit_conversion,
self._cutoff_frequency,
)
@ -180,7 +180,7 @@ class CollisionMatrix(ImagSelfEnergy):
self._triplets_at_q,
self._triplets_map_at_q,
self._ir_map_at_q,
self._temperature,
self._temperature * get_physical_units().KB / get_physical_units().THzToEv,
self._unit_conversion,
self._cutoff_frequency,
)
@ -311,7 +311,11 @@ class CollisionMatrix(ImagSelfEnergy):
freqs = self._frequencies[gp]
sinh = np.where(
freqs > self._cutoff_frequency,
np.sinh(freqs * THzToEv / (2 * Kb * self._temperature)),
np.sinh(
freqs
* get_physical_units().THzToEv
/ (2 * get_physical_units().KB * self._temperature)
),
-1.0,
)
inv_sinh = np.where(sinh > 0, 1.0 / sinh, 0)

9
phono3py/phonon3/gruneisen.py
View File

@ -38,10 +38,10 @@ import sys
import numpy as np
from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import Primitive
from phonopy.structure.grid_points import get_qpoints
from phonopy.units import VaspToTHz
def run_gruneisen_parameters(
@ -133,7 +133,7 @@ class Gruneisen:
nac_params=None,
nac_q_direction=None,
ion_clamped=False,
factor=VaspToTHz,
factor=None,
symprec=1e-5,
):
"""Init method."""
@ -142,7 +142,10 @@ class Gruneisen:
self._scell = supercell
self._pcell = primitive
self._ion_clamped = ion_clamped
self._factor = factor
if factor is None:
self._factor = get_physical_units().DefaultToTHz
else:
self._factor = factor
self._symprec = symprec
self._dm = get_dynamical_matrix(
self._fc2,

22
phono3py/phonon3/imag_self_energy.py
View File

@ -40,7 +40,7 @@ from typing import List, Optional
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.units import EV, Hbar, THz
from phonopy.physical_units import get_physical_units
from phono3py.file_IO import (
write_gamma_detail_to_hdf5,
@ -97,7 +97,11 @@ class ImagSelfEnergy:
# Unit to THz of Gamma
self._unit_conversion = (
18 * np.pi / (Hbar * EV) ** 2 / (2 * np.pi * THz) ** 2 * EV**2
18
* np.pi
/ (get_physical_units().Hbar * get_physical_units().EV) ** 2
/ (2 * np.pi * get_physical_units().THz) ** 2
* get_physical_units().EV ** 2
)
def run(self):
@ -288,7 +292,7 @@ class ImagSelfEnergy:
self._temperature = float(temperature)
def set_temperature(self, temperature):
"""Set temperatures where calculation will be performed."""
"""Set temperature where calculation will be performed."""
warnings.warn(
"Use attribute, ImagSelfEnergy.temperature "
"instead of ImagSelfEnergy.set_temperature().",
@ -387,7 +391,7 @@ class ImagSelfEnergy:
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
self._temperature,
self._temperature * get_physical_units().KB / get_physical_units().THzToEv,
self._g,
_g_zero,
self._cutoff_frequency,
@ -412,7 +416,7 @@ class ImagSelfEnergy:
self._weights_at_q,
self._pp.bz_grid.addresses,
self._frequencies,
self._temperature,
self._temperature * get_physical_units().KB / get_physical_units().THzToEv,
self._g,
_g_zero,
self._cutoff_frequency,
@ -436,7 +440,9 @@ class ImagSelfEnergy:
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
self._temperature,
self._temperature
* get_physical_units().KB
/ get_physical_units().THzToEv,
self._g,
self._g_zero,
self._cutoff_frequency,
@ -471,7 +477,9 @@ class ImagSelfEnergy:
self._weights_at_q,
self._pp.bz_grid.addresses,
self._frequencies,
self._temperature,
self._temperature
* get_physical_units().KB
/ get_physical_units().THzToEv,
g,
_g_zero,
self._cutoff_frequency,

21
phono3py/phonon3/interaction.py
View File

@ -41,9 +41,9 @@ from typing import Literal, Optional, Union
import numpy as np
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
from phonopy.structure.symmetry import Symmetry
from phonopy.units import AMU, EV, Angstrom, Hbar, THz, VaspToTHz
from phono3py.phonon.grid import (
BZGrid,
@ -97,7 +97,7 @@ class Interaction:
fc3: Optional[np.ndarray] = None,
band_indices: Optional[Union[np.ndarray, Sequence]] = None,
constant_averaged_interaction: Optional[float] = None,
frequency_factor_to_THz: float = VaspToTHz,
frequency_factor_to_THz: Optional[float] = None,
frequency_scale_factor: Optional[float] = None,
unit_conversion: Optional[float] = None,
is_mesh_symmetry: bool = True,
@ -115,7 +115,10 @@ class Interaction:
self._band_indices = None
self._set_band_indices(band_indices)
self._constant_averaged_interaction = constant_averaged_interaction
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().DefaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
if fc3 is not None:
@ -125,15 +128,15 @@ class Interaction:
if unit_conversion is None:
num_grid = np.prod(self.mesh_numbers)
self._unit_conversion = (
(Hbar * EV) ** 3
(get_physical_units().Hbar * get_physical_units().EV) ** 3
/ 36
/ 8
* EV**2
/ Angstrom**6
/ (2 * np.pi * THz) ** 3
/ AMU**3
* get_physical_units().EV ** 2
/ get_physical_units().Angstrom ** 6
/ (2 * np.pi * get_physical_units().THz) ** 3
/ get_physical_units().AMU ** 3
/ num_grid
/ EV**2
/ get_physical_units().EV ** 2
)
else:
self._unit_conversion = unit_conversion

13
phono3py/phonon3/joint_dos.py
View File

@ -39,8 +39,8 @@ import warnings
import numpy as np
from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, get_dynamical_matrix
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import Primitive, Supercell
from phonopy.units import VaspToTHz
from phono3py.phonon.func import bose_einstein
from phono3py.phonon.grid import BZGrid, get_grid_point_from_address
@ -66,7 +66,7 @@ class JointDos:
sigma=None,
sigma_cutoff=None,
cutoff_frequency=None,
frequency_factor_to_THz=VaspToTHz,
frequency_factor_to_THz=None,
frequency_scale_factor=1.0,
is_mesh_symmetry=True,
symprec=1e-5,
@ -90,7 +90,10 @@ class JointDos:
self._cutoff_frequency = 0
else:
self._cutoff_frequency = cutoff_frequency
self._frequency_factor_to_THz = frequency_factor_to_THz
if frequency_factor_to_THz is None:
self._frequency_factor_to_THz = get_physical_units().DefaultToTHz
else:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
self._is_mesh_symmetry = is_mesh_symmetry
self._symprec = symprec
@ -388,7 +391,9 @@ class JointDos:
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
self._temperature,
self._temperature
* get_physical_units().KB
/ get_physical_units().THzToEv,
g,
self._g_zero,
self._cutoff_frequency,

15
phono3py/phonon3/real_self_energy.py
View File

@ -38,7 +38,7 @@ import sys
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.units import EV, Hbar, THz
from phonopy.physical_units import get_physical_units
from phono3py.file_IO import (
write_real_self_energy_at_grid_point,
@ -115,7 +115,12 @@ class RealSelfEnergy:
self._real_self_energies = None
# Unit to THz of Delta
self._unit_conversion = 18 / (Hbar * EV) ** 2 / (2 * np.pi * THz) ** 2 * EV**2
self._unit_conversion = (
18
/ (get_physical_units().Hbar * get_physical_units().EV) ** 2
/ (2 * np.pi * get_physical_units().THz) ** 2
* get_physical_units().EV ** 2
)
def run(self):
"""Calculate real-part of self-energies."""
@ -240,7 +245,7 @@ class RealSelfEnergy:
self._weights_at_q,
self._frequencies,
self._band_indices,
self._temperature,
self._temperature * get_physical_units().KB / get_physical_units().THzToEv,
self._epsilon,
self._unit_conversion,
self._cutoff_frequency,
@ -277,7 +282,9 @@ class RealSelfEnergy:
self._weights_at_q,
self._frequencies,
self._band_indices,
self._temperature,
self._temperature
* get_physical_units().KB
/ get_physical_units().THzToEv,
self._epsilon,
self._unit_conversion,
self._cutoff_frequency,

11
phono3py/sscha/sscha.py
View File

@ -43,7 +43,7 @@ Formulae implemented are based on these papers:
import numpy as np
from phonopy.harmonic.dynmat_to_fc import DynmatToForceConstants
from phonopy.units import VaspToTHz
from phonopy.physical_units import get_physical_units
from phono3py.phonon.func import sigma_squared
@ -97,7 +97,7 @@ class SupercellPhonon:
"""
def __init__(self, supercell, force_constants, frequency_factor_to_THz=VaspToTHz):
def __init__(self, supercell, force_constants, frequency_factor_to_THz=None):
"""Init method.
Parameters
@ -107,7 +107,7 @@ class SupercellPhonon:
force_constants : array_like
Second order force constants.
shape=(num_satom, num_satom, 3, 3), dtype='double', order='C'
frequency_factor_to_THz : float
frequency_factor_to_THz : float, optional
Frequency conversion factor to THz.
"""
@ -124,7 +124,10 @@ class SupercellPhonon:
)
eigvals, eigvecs = np.linalg.eigh(dynmat)
freqs = np.sqrt(np.abs(eigvals)) * np.sign(eigvals)
freqs *= frequency_factor_to_THz
if frequency_factor_to_THz is None:
freqs *= get_physical_units().DefaultToTHz
else:
freqs *= frequency_factor_to_THz
self._eigenvalues = np.array(eigvals, dtype="double", order="C")
self._eigenvectors = np.array(eigvecs, dtype="double", order="C")
self._frequencies = np.array(freqs, dtype="double", order="C")

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.0"
__version__ = "3.15.1"

2
pyproject.toml
View File

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

6
test/phonon/test_velocity_operator.py
View File

@ -2,7 +2,7 @@
import numpy as np
from phonopy import Phonopy
from phonopy.units import THzToCm, VaspToTHz
from phonopy.physical_units import get_physical_units
from phono3py.phonon.velocity_operator import VelocityOperator
@ -194,7 +194,7 @@ def test_gv_operator_nacl(ph_nacl: Phonopy):
eigvals, eigvecs = np.linalg.eigh(dm)
np.testing.assert_allclose(
eigvals * VaspToTHz * THzToCm,
eigvals * get_physical_units().DefaultToTHz * get_physical_units().THzToCm,
eigvals_NaCl_Ref,
atol=0.00001,
rtol=0.00001,
@ -406,7 +406,7 @@ def test_gv_operator_si(ph_si: Phonopy):
eigvals, eigvecs = np.linalg.eigh(dm)
np.testing.assert_allclose(
eigvals * VaspToTHz * THzToCm,
eigvals * get_physical_units().DefaultToTHz * get_physical_units().THzToCm,
eigvals_si_Ref,
atol=0.00001,
rtol=0.00001,