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

Merge branch 'develop' of github.com:phonopy/phono3py into develop

This commit is contained in:
Atsushi Togo 2025-06-17 15:59:35 +09:00
parent 518b6feebf 685d5f5e55
commit 36ae0c1ffb
4 changed files with 77 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
c/interaction.c
View File

@ -56,9 +56,10 @@ static void real_to_normal(
const lapack_complex_double *eigvecs2, const double *fc3,
const int64_t is_compact_fc3, const double q_vecs[3][3], /* q0, q1, q2 */
const AtomTriplets *atom_triplets, const double *masses,
const int64_t *band_indices, const int64_t num_band,
const double cutoff_frequency, const int64_t triplet_index,
const int64_t num_triplets, const int64_t openmp_per_triplets);
const int64_t *band_indices, const int64_t num_band0,
const int64_t num_band, const double cutoff_frequency,
const int64_t triplet_index, const int64_t num_triplets,
const int64_t openmp_per_triplets);
static void real_to_normal_sym_q(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, double *const freqs[3],
@ -79,7 +80,7 @@ void itr_get_interaction(
const AtomTriplets *atom_triplets, const double *masses,
const int64_t *band_indices, const int64_t symmetrize_fc3_q,
const double cutoff_frequency, const int64_t openmp_per_triplets) {
int64_t(*g_pos)[4];
int64_t (*g_pos)[4];
int64_t i;
int64_t num_band, num_band0, num_band_prod, num_g_pos;
@ -90,11 +91,11 @@ void itr_get_interaction(
num_band_prod = num_band0 * num_band * num_band;
#ifdef _OPENMP
#pragma omp parallel for schedule(guided) private( \
num_g_pos, g_pos) if (openmp_per_triplets)
#pragma omp parallel for schedule(guided) \
private(num_g_pos, g_pos) if (openmp_per_triplets)
#endif
for (i = 0; i < num_triplets; i++) {
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);
num_g_pos = ise_set_g_pos(g_pos, num_band0, num_band,
g_zero + i * num_band_prod);
@ -177,7 +178,7 @@ void itr_get_interaction_at_triplet(
eigenvectors + triplet[1] * num_band * num_band,
eigenvectors + triplet[2] * num_band * num_band, fc3,
is_compact_fc3, q_vecs, /* q0, q1, q2 */
atom_triplets, masses, band_indices, num_band,
atom_triplets, masses, band_indices, num_band0, num_band,
cutoff_frequency, triplet_index, num_triplets,
openmp_per_triplets);
}
@ -191,9 +192,10 @@ static void real_to_normal(
const lapack_complex_double *eigvecs2, const double *fc3,
const int64_t is_compact_fc3, const double q_vecs[3][3], /* q0, q1, q2 */
const AtomTriplets *atom_triplets, const double *masses,
const int64_t *band_indices, const int64_t num_band,
const double cutoff_frequency, const int64_t triplet_index,
const int64_t num_triplets, const int64_t openmp_per_triplets) {
const int64_t *band_indices, const int64_t num_band0,
const int64_t num_band, const double cutoff_frequency,
const int64_t triplet_index, const int64_t num_triplets,
const int64_t openmp_per_triplets) {
lapack_complex_double *fc3_reciprocal;
lapack_complex_double comp_zero;
int64_t i;
@ -215,8 +217,8 @@ static void real_to_normal(
#endif
reciprocal_to_normal_squared(
fc3_normal_squared, g_pos, num_g_pos, fc3_reciprocal, freqs0, freqs1,
freqs2, eigvecs0, eigvecs1, eigvecs2, masses, band_indices, num_band,
cutoff_frequency, openmp_per_triplets);
freqs2, eigvecs0, eigvecs1, eigvecs2, masses, band_indices, num_band0,
num_band, cutoff_frequency, openmp_per_triplets);
free(fc3_reciprocal);
fc3_reciprocal = NULL;
@ -256,8 +258,8 @@ static void real_to_normal_sym_q(
freqs[index_exchange[i][2]], eigvecs[index_exchange[i][0]],
eigvecs[index_exchange[i][1]], eigvecs[index_exchange[i][2]], fc3,
is_compact_fc3, q_vecs_ex, /* q0, q1, q2 */
atom_triplets, masses, band_indices, num_band, cutoff_frequency,
triplet_index, num_triplets, openmp_per_triplets);
atom_triplets, masses, band_indices, num_band0, num_band,
cutoff_frequency, triplet_index, num_triplets, openmp_per_triplets);
for (j = 0; j < num_band0; j++) {
for (k = 0; k < num_band; k++) {
for (l = 0; l < num_band; l++) {

138
c/reciprocal_to_normal.c
View File

@ -49,21 +49,18 @@
#include "lapack_wrapper.h"
static void reciprocal_to_normal_squared_no_threading(
static void get_fc3_e0_e1_e2(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *e0, const lapack_complex_double *e1,
const lapack_complex_double *e2, const int64_t *band_indices,
const int64_t num_band, const double cutoff_frequency);
static void reciprocal_to_normal_squared_g_threading(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *e0, const lapack_complex_double *e1,
const lapack_complex_double *e2, const int64_t *band_indices,
const int64_t num_band, const double cutoff_frequency);
const int64_t num_band0, const int64_t num_band,
const double cutoff_frequency, const int64_t openmp_per_triplets);
static void get_fc3_e0(lapack_complex_double *fc3_e0,
const lapack_complex_double *fc3_reciprocal,
const lapack_complex_double *e0,
const int64_t band_index_0, const int64_t num_band);
static double get_fc3_sum(const lapack_complex_double *e1,
const lapack_complex_double *e2,
const lapack_complex_double *fc3_reciprocal,
@ -89,8 +86,9 @@ void reciprocal_to_normal_squared(
const lapack_complex_double *eigvecs0,
const lapack_complex_double *eigvecs1,
const lapack_complex_double *eigvecs2, const double *masses,
const int64_t *band_indices, const int64_t num_band,
const double cutoff_frequency, const int64_t openmp_per_triplets) {
const int64_t *band_indices, const int64_t num_band0,
const int64_t num_band, const double cutoff_frequency,
const int64_t openmp_per_triplets) {
int64_t i, j, ij, num_atom;
double *inv_sqrt_masses;
lapack_complex_double *e0, *e1, *e2;
@ -144,17 +142,10 @@ void reciprocal_to_normal_squared(
free(inv_sqrt_masses);
inv_sqrt_masses = NULL;
if (openmp_per_triplets) {
reciprocal_to_normal_squared_no_threading(
fc3_normal_squared, g_pos, num_g_pos, fc3_reciprocal, freqs0,
freqs1, freqs2, e0, e1, e2, band_indices, num_band,
cutoff_frequency);
} else {
reciprocal_to_normal_squared_g_threading(
fc3_normal_squared, g_pos, num_g_pos, fc3_reciprocal, freqs0,
freqs1, freqs2, e0, e1, e2, band_indices, num_band,
cutoff_frequency);
}
get_fc3_e0_e1_e2(fc3_normal_squared, g_pos, num_g_pos, fc3_reciprocal,
freqs0, freqs1, freqs2, e0, e1, e2, band_indices,
num_band0, num_band, cutoff_frequency,
openmp_per_triplets);
free(e0);
e0 = NULL;
@ -162,54 +153,45 @@ void reciprocal_to_normal_squared(
e2 = NULL;
}
// This function is more efficient than the one with multithreading
// getting rid of no concurrency over g.
static void reciprocal_to_normal_squared_no_threading(
// This is less efficient than the one without multithreading
// but can be called when unit cell is large.
static void get_fc3_e0_e1_e2(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *e0, const lapack_complex_double *e1,
const lapack_complex_double *e2, const int64_t *band_indices,
const int64_t num_band, const double cutoff_frequency) {
int64_t i, j, k, ll, bi_prev;
lapack_complex_double *fc3_e0, fc3_elem, zero;
const int64_t num_band0, const int64_t num_band,
const double cutoff_frequency, const int64_t openmp_per_triplets) {
int64_t i;
lapack_complex_double *fc3_e0, zero;
zero = lapack_make_complex_double(0, 0);
bi_prev = -1;
fc3_e0 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) *
num_band * num_band);
num_band0 * num_band * num_band);
for (i = 0; i < num_band0 * num_band * num_band; i++) {
fc3_e0[i] = zero;
}
#ifdef _OPENMP
#pragma omp parallel for if (!openmp_per_triplets)
#endif
for (i = 0; i < num_band0; i++) {
get_fc3_e0(fc3_e0 + i * num_band * num_band, fc3_reciprocal, e0,
band_indices[i], num_band);
}
#ifdef _OPENMP
#pragma omp parallel for if (!openmp_per_triplets)
#endif
for (i = 0; i < num_g_pos; i++) {
if (freqs0[band_indices[g_pos[i][0]]] > cutoff_frequency &&
freqs1[g_pos[i][1]] > cutoff_frequency &&
freqs2[g_pos[i][2]] > cutoff_frequency) {
if (bi_prev != g_pos[i][0]) {
bi_prev = g_pos[i][0];
for (j = 0; j < num_band * num_band; j++) {
fc3_e0[j] = zero;
}
for (j = 0; j < num_band; j++) {
for (k = 0; k < num_band; k++) {
for (ll = 0; ll < num_band; ll++) {
fc3_elem = phonoc_complex_prod(
fc3_reciprocal[j * num_band * num_band +
k * num_band + ll],
e0[band_indices[g_pos[i][0]] * num_band + j]);
fc3_e0[k * num_band + ll] =
lapack_make_complex_double(
lapack_complex_double_real(
fc3_e0[k * num_band + ll]) +
lapack_complex_double_real(fc3_elem),
lapack_complex_double_imag(
fc3_e0[k * num_band + ll]) +
lapack_complex_double_imag(fc3_elem));
}
}
}
}
fc3_normal_squared[g_pos[i][3]] =
get_fc3_sum(e1 + g_pos[i][1] * num_band,
e2 + g_pos[i][2] * num_band, fc3_e0, num_band) /
get_fc3_sum(
e1 + g_pos[i][1] * num_band, e2 + g_pos[i][2] * num_band,
fc3_e0 + g_pos[i][0] * num_band * num_band, num_band) /
(freqs0[band_indices[g_pos[i][0]]] * freqs1[g_pos[i][1]] *
freqs2[g_pos[i][2]]);
} else {
@ -221,33 +203,23 @@ static void reciprocal_to_normal_squared_no_threading(
fc3_e0 = NULL;
}
// This is less efficient than the one without multithreading
// but can be called when unit cell is large.
static void reciprocal_to_normal_squared_g_threading(
double *fc3_normal_squared, const int64_t (*g_pos)[4],
const int64_t num_g_pos, const lapack_complex_double *fc3_reciprocal,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *e0, const lapack_complex_double *e1,
const lapack_complex_double *e2, const int64_t *band_indices,
const int64_t num_band, const double cutoff_frequency) {
int64_t i;
static void get_fc3_e0(lapack_complex_double *fc3_e0,
const lapack_complex_double *fc3_reciprocal,
const lapack_complex_double *e0,
const int64_t band_index_0, const int64_t num_band) {
int64_t j, k;
lapack_complex_double fc3_elem;
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_g_pos; i++) {
if (freqs0[band_indices[g_pos[i][0]]] > cutoff_frequency &&
freqs1[g_pos[i][1]] > cutoff_frequency &&
freqs2[g_pos[i][2]] > cutoff_frequency) {
fc3_normal_squared[g_pos[i][3]] =
get_fc3_sum_blas_like(e0 + band_indices[g_pos[i][0]] * num_band,
e1 + g_pos[i][1] * num_band,
e2 + g_pos[i][2] * num_band,
fc3_reciprocal, num_band) /
(freqs0[band_indices[g_pos[i][0]]] * freqs1[g_pos[i][1]] *
freqs2[g_pos[i][2]]);
} else {
fc3_normal_squared[g_pos[i][3]] = 0;
for (j = 0; j < num_band; j++) {
for (k = 0; k < num_band * num_band; k++) {
fc3_elem =
phonoc_complex_prod(fc3_reciprocal[j * num_band * num_band + k],
e0[band_index_0 * num_band + j]);
fc3_e0[k] = lapack_make_complex_double(
lapack_complex_double_real(fc3_e0[k]) +
lapack_complex_double_real(fc3_elem),
lapack_complex_double_imag(fc3_e0[k]) +
lapack_complex_double_imag(fc3_elem));
}
}
}

5
c/reciprocal_to_normal.h
View File

@ -46,7 +46,8 @@ void reciprocal_to_normal_squared(
const lapack_complex_double *eigvecs0,
const lapack_complex_double *eigvecs1,
const lapack_complex_double *eigvecs2, const double *masses,
const int64_t *band_indices, const int64_t num_band,
const double cutoff_frequency, const int64_t openmp_per_triplets);
const int64_t *band_indices, const int64_t num_band0,
const int64_t num_band, const double cutoff_frequency,
const int64_t openmp_per_triplets);
#endif

4
phono3py/conductivity/rta_base.py
View File

@ -204,7 +204,7 @@ class ConductivityRTABase(ConductivityBase):
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)
print("Number of triplets: %d" % num_triplets, flush=True)
if (
self._is_full_pp
@ -288,7 +288,7 @@ class ConductivityRTABase(ConductivityBase):
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),
((num_temp,) + self._pp.interaction_strength.shape),
dtype="double",
order="C",
)