Merge pull request #199 from phonopy/r0-ave-with-atom-triplets

fc3 transform with r0-average and all-shortest atom triplets

c/_phono3py.c

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

c/phono3py.c

@@ -67,8 +67,8 @@ long ph3py_get_interaction(
     const long multi_dims[2], const long (*multiplicity)[2],
     const double *masses, const long *p2s_map, const long *s2p_map,
     const long *band_indices, const long symmetrize_fc3_q,
-    const long make_r0_average, const double cutoff_frequency,
-    const long openmp_per_triplets) {
+    const long make_r0_average, const char *all_shortest,
+    const double cutoff_frequency, const long openmp_per_triplets) {
     ConstBZGrid *bzgrid;
     AtomTriplets *atom_triplets;
     long i, j;
@@ -100,6 +100,7 @@ long ph3py_get_interaction(
     atom_triplets->p2s_map = p2s_map;
     atom_triplets->s2p_map = s2p_map;
     atom_triplets->make_r0_average = make_r0_average;
+    atom_triplets->all_shortest = all_shortest;
 
     itr_get_interaction(fc3_normal_squared, g_zero, frequencies,
                         (lapack_complex_double *)eigenvectors, triplets,
@@ -129,7 +130,8 @@ long ph3py_get_pp_collision(
     const double *masses, const long *p2s_map, const long *s2p_map,
     const Larray *band_indices, const Darray *temperatures, const long is_NU,
     const long symmetrize_fc3_q, const long make_r0_average,
-    const double cutoff_frequency, const long openmp_per_triplets) {
+    const char *all_shortest, const double cutoff_frequency,
+    const long openmp_per_triplets) {
     ConstBZGrid *bzgrid;
     AtomTriplets *atom_triplets;
     long i, j;
@@ -163,6 +165,7 @@ long ph3py_get_pp_collision(
     atom_triplets->p2s_map = p2s_map;
     atom_triplets->s2p_map = s2p_map;
     atom_triplets->make_r0_average = make_r0_average;
+    atom_triplets->all_shortest = all_shortest;
 
     ppc_get_pp_collision(imag_self_energy, relative_grid_address, frequencies,
                          (lapack_complex_double *)eigenvectors, triplets,
@@ -191,7 +194,8 @@ long ph3py_get_pp_collision_with_sigma(
     const double *masses, const long *p2s_map, const long *s2p_map,
     const Larray *band_indices, const Darray *temperatures, const long is_NU,
     const long symmetrize_fc3_q, const long make_r0_average,
-    const double cutoff_frequency, const long openmp_per_triplets) {
+    const char *all_shortest, const double cutoff_frequency,
+    const long openmp_per_triplets) {
     ConstBZGrid *bzgrid;
     AtomTriplets *atom_triplets;
     long i, j;
@@ -223,6 +227,7 @@ long ph3py_get_pp_collision_with_sigma(
     atom_triplets->p2s_map = p2s_map;
     atom_triplets->s2p_map = s2p_map;
     atom_triplets->make_r0_average = make_r0_average;
+    atom_triplets->all_shortest = all_shortest;
 
     ppc_get_pp_collision_with_sigma(
         imag_self_energy, sigma, sigma_cutoff, frequencies,

c/phono3py.h

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

c/real_to_reciprocal.c

@@ -54,7 +54,7 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
                                         const long is_compact_fc3,
                                         const AtomTriplets *atom_triplets,
                                         const long pi0, const long pi1,
-                                        const long pi2);
+                                        const long pi2, const long leg_index);
 static lapack_complex_double get_phase_factor(const double q[3],
                                               const double (*svecs)[3],
                                               const long multi[2]);
@@ -100,10 +100,12 @@ static void real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
         i = ijk / (num_patom * num_patom);
         j = (ijk - (i * num_patom * num_patom)) / num_patom;
         k = ijk % num_patom;
+
         pre_phase_factor = get_pre_phase_factor(i, q_vecs, atom_triplets);
 
         real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[2], fc3,
-                                    is_compact_fc3, atom_triplets, i, j, k);
+                                    is_compact_fc3, atom_triplets, i, j, k,
+                                    atom_triplets->make_r0_average);
         for (l = 0; l < 3; l++) {
             for (m = 0; m < 3; m++) {
                 for (n = 0; n < 3; n++) {
@@ -121,7 +123,8 @@ static void real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
 
         if (atom_triplets->make_r0_average) {
             real_to_reciprocal_elements(fc3_rec_elem, q_vecs[0], q_vecs[2], fc3,
-                                        is_compact_fc3, atom_triplets, i, j, k);
+                                        is_compact_fc3, atom_triplets, i, j, k,
+                                        2);
             for (l = 0; l < 3; l++) {
                 for (m = 0; m < 3; m++) {
                     for (n = 0; n < 3; n++) {
@@ -141,7 +144,8 @@ static void real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
             }
 
             real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[0], fc3,
-                                        is_compact_fc3, atom_triplets, i, j, k);
+                                        is_compact_fc3, atom_triplets, i, j, k,
+                                        3);
             for (l = 0; l < 3; l++) {
                 for (m = 0; m < 3; m++) {
                     for (n = 0; n < 3; n++) {
@@ -169,7 +173,7 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
                                         const long is_compact_fc3,
                                         const AtomTriplets *atom_triplets,
                                         const long pi0, const long pi1,
-                                        const long pi2) {
+                                        const long pi2, const long leg_index) {
     long i, j, k, l;
     long num_satom, adrs_shift, adrs_vec1, adrs_vec2;
     lapack_complex_double phase_factor, phase_factor1, phase_factor2;
@@ -200,18 +204,40 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
             if (atom_triplets->s2p_map[k] != atom_triplets->p2s_map[pi2]) {
                 continue;
             }
+            if (leg_index > 1) {
+                if (atom_triplets->all_shortest[pi0 * num_satom * num_satom +
+                                                j * num_satom + k]) {
+                    continue;
+                }
+            }
             adrs_vec2 = k * atom_triplets->multi_dims[1] + pi0;
             phase_factor2 =
                 get_phase_factor(q2, atom_triplets->svecs,
                                  atom_triplets->multiplicity[adrs_vec2]);
             adrs_shift =
                 i * 27 * num_satom * num_satom + j * 27 * num_satom + k * 27;
             phase_factor = phonoc_complex_prod(phase_factor1, phase_factor2);
-            for (l = 0; l < 27; l++) {
-                fc3_rec_real[l] += lapack_complex_double_real(phase_factor) *
-                                   fc3[adrs_shift + l];
-                fc3_rec_imag[l] += lapack_complex_double_imag(phase_factor) *
-                                   fc3[adrs_shift + l];
+
+            if ((leg_index == 1) &&
+                (atom_triplets->all_shortest[pi0 * num_satom * num_satom +
+                                             j * num_satom + k])) {
+                for (l = 0; l < 27; l++) {
+                    fc3_rec_real[l] +=
+                        lapack_complex_double_real(phase_factor) *
+                        fc3[adrs_shift + l] * 3;
+                    fc3_rec_imag[l] +=
+                        lapack_complex_double_imag(phase_factor) *
+                        fc3[adrs_shift + l] * 3;
+                }
+            } else {
+                for (l = 0; l < 27; l++) {
+                    fc3_rec_real[l] +=
+                        lapack_complex_double_real(phase_factor) *
+                        fc3[adrs_shift + l];
+                    fc3_rec_imag[l] +=
+                        lapack_complex_double_imag(phase_factor) *
+                        fc3[adrs_shift + l];
+                }
             }
         }
     }
@@ -227,7 +253,7 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
 static lapack_complex_double get_pre_phase_factor(
     const long i_patom, const double q_vecs[3][3],
     const AtomTriplets *atom_triplets) {
-    long i, j, svecs_adrs;
+    long j, svecs_adrs;
     double pre_phase;
     lapack_complex_double pre_phase_factor;
 

c/real_to_reciprocal.h

@@ -45,6 +45,7 @@ typedef struct {
     const long *p2s_map;
     const long *s2p_map;
     long make_r0_average;
+    const char *all_shortest;
 } AtomTriplets;
 
 void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,

phono3py/conductivity/direct_solution.py

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