From 1d0e2a55ddfe5e41e4e8f6e135e19a5cdeddacf3 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Thu, 10 Mar 2016 10:57:45 +0900 Subject: [PATCH] OpenMP support for the code from real-space fc3 to reciprocal-space fc3 --- anharmonic/cui/phono3py_argparse.py | 20 +++-- anharmonic/cui/settings.py | 19 +++++ anharmonic/cui/show_log.py | 2 +- anharmonic/phonon3/__init__.py | 2 + anharmonic/phonon3/conductivity_RTA.py | 2 + c/anharmonic/phonon3/interaction.c | 3 +- c/anharmonic/phonon3/real_to_reciprocal.c | 29 ++++++- c/anharmonic/phonon3/reciprocal_to_normal.c | 85 +++++++++---------- c/anharmonic_h/phonon3_h/real_to_reciprocal.h | 3 +- scripts/phono3py | 1 + 10 files changed, 107 insertions(+), 59 deletions(-) diff --git a/anharmonic/cui/phono3py_argparse.py b/anharmonic/cui/phono3py_argparse.py index 25d347b6..39d8816f 100644 --- a/anharmonic/cui/phono3py_argparse.py +++ b/anharmonic/cui/phono3py_argparse.py @@ -36,7 +36,8 @@ from optparse import OptionParser def get_parser(): parser = OptionParser() - parser.set_defaults(band_indices=None, + parser.set_defaults(average_pp_interaction=False, + band_indices=None, band_paths=None, band_points=None, cell_filename=None, @@ -91,9 +92,13 @@ def get_parser(): mesh_numbers=None, mesh_divisors=None, no_kappa_stars=False, + output_filename=None, phonon_supercell_dimension=None, pinv_cutoff=1.0e-8, + pp_unit_conversion=None, primitive_axis=None, + qpoints=None, + quiet=False, q_direction=None, read_amplitude=False, read_collision=None, @@ -101,10 +106,6 @@ def get_parser(): read_fc3=False, read_gamma=False, run_with_g=True, - output_filename=None, - pp_unit_conversion=None, - qpoints=None, - quiet=False, scattering_event_class=None, show_num_triplets=False, sigma=None, @@ -116,8 +117,8 @@ def get_parser(): tstep=None, tsym_type=None, uplo='L', - average_pp_interaction=False, verbose=False, + with_g_zero=False, write_amplitude=False, write_collision=False, write_detailed_gamma=False, @@ -385,6 +386,13 @@ def get_parser(): parser.add_option( "-v", "--verbose", dest="verbose", action="store_true", help="Detailed run-time information is displayed") + parser.add_option( + "--with_g_zero", dest="with_g_zero", + action="store_true", + help=("Avoid calculating ph-ph interaction for RTA conductivity when " + "integration weight is zero. This is activated when only using " + "tetrahedron method and no smearing method is used. Average " + "ph-ph interaction is not obtained with this option.")) parser.add_option( "--wgp", "--write_grid_points", dest="write_grid_points", action="store_true", diff --git a/anharmonic/cui/settings.py b/anharmonic/cui/settings.py index 3e396af1..5d0e3914 100644 --- a/anharmonic/cui/settings.py +++ b/anharmonic/cui/settings.py @@ -46,6 +46,7 @@ class Phono3pySettings(Settings): self._pp_conversion_factor = None self._scattering_event_class = None # scattering event class 1 or 2 self._temperatures = None + self._with_g_zero = False self._write_amplitude = False self._write_collision = False self._write_detailed_gamma = False @@ -292,6 +293,12 @@ class Phono3pySettings(Settings): def get_average_pp_interaction(self): return self._average_pp_interaction + def set_with_g_zero(self, with_g_zero): + self._with_g_zero = with_g_zero + + def get_with_g_zero(self): + return self._with_g_zero + def set_write_amplitude(self, write_amplitude): self._write_amplitude = write_amplitude @@ -488,6 +495,10 @@ class Phono3pyConfParser(ConfParser): if self._options.average_pp_interaction: self._confs['average_pp_interaction'] = '.true.' + if opt.dest == 'with_g_zero': + if self._options.with_g_zero: + self._confs['with_g_zero'] = '.true.' + if opt.dest == 'write_amplitude': if self._options.write_amplitude: self._confs['write_amplitude'] = '.true.' @@ -710,6 +721,10 @@ class Phono3pyConfParser(ConfParser): if confs['average_pp_interaction'] == '.true.': self.set_parameter('average_pp_interaction', True) + if conf_key == 'with_g_zero': + if confs['with_g_zero'] == '.true.': + self.set_parameter('with_g_zero', True) + if conf_key == 'write_amplitude': if confs['write_amplitude'] == '.true.': self.set_parameter('write_amplitude', True) @@ -905,6 +920,10 @@ class Phono3pyConfParser(ConfParser): self._settings.set_average_pp_interaction( params['average_pp_interaction']) + # Avoid calculating useless ph-ph calculation for RTA conductivity + if params.has_key('with_g_zero'): + self._settings.set_with_g_zero(params['with_g_zero']) + # Write phonon-phonon interaction amplitudes to hdf5 if params.has_key('write_amplitude'): self._settings.set_write_amplitude(params['write_amplitude']) diff --git a/anharmonic/cui/show_log.py b/anharmonic/cui/show_log.py index c02e78f1..d986cb9a 100644 --- a/anharmonic/cui/show_log.py +++ b/anharmonic/cui/show_log.py @@ -203,7 +203,7 @@ def show_phono3py_settings(settings, if (settings.get_average_pp_interaction() and (settings.get_is_bterta() or settings.get_is_lbte())): print("Use averaged ph-ph interaction") - + if log_level > 1: print("Frequency factor to THz: %s" % frequency_factor_to_THz) if frequency_step is not None: diff --git a/anharmonic/phonon3/__init__.py b/anharmonic/phonon3/__init__.py index 85c067f7..9a361fd4 100644 --- a/anharmonic/phonon3/__init__.py +++ b/anharmonic/phonon3/__init__.py @@ -428,6 +428,7 @@ class Phono3py: is_kappa_star=True, gv_delta_q=None, # for group velocity run_with_g=True, # integration weights for smearing method, too + with_g_zero=False, pinv_cutoff=1.0e-8, # for pseudo-inversion of collision matrix write_gamma=False, read_gamma=False, @@ -476,6 +477,7 @@ class Phono3py: is_kappa_star=is_kappa_star, gv_delta_q=gv_delta_q, run_with_g=run_with_g, + with_g_zero=with_g_zero, write_gamma=write_gamma, read_gamma=read_gamma, input_filename=input_filename, diff --git a/anharmonic/phonon3/conductivity_RTA.py b/anharmonic/phonon3/conductivity_RTA.py index d330fb92..cc52329b 100644 --- a/anharmonic/phonon3/conductivity_RTA.py +++ b/anharmonic/phonon3/conductivity_RTA.py @@ -24,6 +24,7 @@ def get_thermal_conductivity_RTA( is_kappa_star=True, gv_delta_q=1e-4, run_with_g=True, # integration weights from gaussian smearing function + with_g_zero=False, write_gamma=False, read_gamma=False, input_filename=None, @@ -52,6 +53,7 @@ def get_thermal_conductivity_RTA( is_kappa_star=is_kappa_star, gv_delta_q=gv_delta_q, run_with_g=run_with_g, + with_g_zero=with_g_zero, log_level=log_level) if read_gamma: diff --git a/c/anharmonic/phonon3/interaction.c b/c/anharmonic/phonon3/interaction.c index 7ef6e19d..3f2822f8 100644 --- a/c/anharmonic/phonon3/interaction.c +++ b/c/anharmonic/phonon3/interaction.c @@ -300,7 +300,8 @@ static void real_to_normal(double *fc3_normal_squared, shortest_vectors, multiplicity, p2s_map, - s2p_map); + s2p_map, + openmp_at_bands); if (openmp_at_bands) { #ifdef MEASURE_R2N diff --git a/c/anharmonic/phonon3/real_to_reciprocal.c b/c/anharmonic/phonon3/real_to_reciprocal.c index 43501bd2..c638a1a4 100644 --- a/c/anharmonic/phonon3/real_to_reciprocal.c +++ b/c/anharmonic/phonon3/real_to_reciprocal.c @@ -59,17 +59,21 @@ void real_to_reciprocal(lapack_complex_double *fc3_reciprocal, const Darray *shortest_vectors, const Iarray *multiplicity, const int *p2s_map, - const int *s2p_map) + const int *s2p_map, + const int openmp_at_bands) { - int i, j, k, num_patom; + int i, j, k, jk, num_patom; double pre_phase; lapack_complex_double pre_phase_factor; num_patom = multiplicity->dims[1]; for (i = 0; i < num_patom; i++) { - for (j = 0; j < num_patom; j++) { - for (k = 0; k < num_patom; k++) { + if (openmp_at_bands) { +#pragma omp parallel for private(j, k) + for (jk = 0; jk < num_patom * num_patom; jk++) { + j = jk / num_patom; + k = jk % num_patom; real_to_reciprocal_elements(fc3_reciprocal + i * 27 * num_patom * num_patom + j * 27 * num_patom + @@ -83,6 +87,23 @@ void real_to_reciprocal(lapack_complex_double *fc3_reciprocal, i, j, k); } + } else { + for (j = 0; j < num_patom; j++) { + for (k = 0; k < num_patom; k++) { + real_to_reciprocal_elements(fc3_reciprocal + + i * 27 * num_patom * num_patom + + j * 27 * num_patom + + k * 27, + q, + fc3, + shortest_vectors, + multiplicity, + p2s_map, + s2p_map, + i, j, k); + + } + } } pre_phase = 0; diff --git a/c/anharmonic/phonon3/reciprocal_to_normal.c b/c/anharmonic/phonon3/reciprocal_to_normal.c index 460cc7c0..89c9c021 100644 --- a/c/anharmonic/phonon3/reciprocal_to_normal.c +++ b/c/anharmonic/phonon3/reciprocal_to_normal.c @@ -32,6 +32,7 @@ /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */ /* POSSIBILITY OF SUCH DAMAGE. */ +#include #include #include #include @@ -134,7 +135,8 @@ void reciprocal_to_normal_squared_openmp const int num_band, const double cutoff_frequency) { - int i, j, k, jk, bi, num_atom; + int i, j, k, l, ijk, jk, bi, num_atom, count; + int *n; #ifdef MEASURE_R2N double loopTotalCPUTime, loopTotalWallTime; @@ -144,61 +146,52 @@ void reciprocal_to_normal_squared_openmp num_atom = num_band / 3; - for (i = 0; i < num_band0; i++) { - bi = band_indices[i]; - if (freqs0[bi] > cutoff_frequency) { + n = (int*)malloc(sizeof(int) * num_band0 * num_band * num_band); + + count = 0; + for (ijk = 0; ijk < num_band0 * num_band * num_band; ijk++) { + if (g_zero[ijk]) { + fc3_normal_squared[ijk] = 0; + } else { + n[count] = ijk; + count++; + } + } #ifdef MEASURE_R2N - loopStartWallTime = time(NULL); - loopStartCPUTime = clock(); + loopStartWallTime = time(NULL); + loopStartCPUTime = clock(); #endif -#pragma omp parallel for private(j, k) - for (jk = 0; jk < num_band * num_band; jk++) { - if (g_zero[i * num_band * num_band + jk]) { - fc3_normal_squared[i * num_band * num_band + jk] = 0; - continue; - } - j = jk / num_band; - k = jk % num_band; - fc3_normal_squared[i * num_band * num_band + jk] = - get_fc3_sum(i, j, k, bi, - freqs0, freqs1, freqs2, - eigvecs0, eigvecs1, eigvecs2, - fc3_reciprocal, - masses, - num_atom, - num_band, - cutoff_frequency); - } +#pragma omp parallel for private(ijk, i, j, k, jk, bi) + for (l = 0; l < count; l++) { + ijk = n[l]; + i = ijk / (num_band * num_band); + jk = ijk % (num_band * num_band); + j = jk / num_band; + k = jk % num_band; + bi = band_indices[i]; + if (freqs0[bi] > cutoff_frequency) { + fc3_normal_squared[n[l]] = get_fc3_sum(i, j, k, bi, + freqs0, freqs1, freqs2, + eigvecs0, eigvecs1, eigvecs2, + fc3_reciprocal, + masses, + num_atom, + num_band, + cutoff_frequency); + } else { + fc3_normal_squared[n[l]] = 0; + } + } #ifdef MEASURE_R2N loopTotalCPUTime = (double)(clock() - loopStartCPUTime) / CLOCKS_PER_SEC; loopTotalWallTime = difftime(time(NULL), loopStartWallTime); - printf(" Band index %d/%d %1.3fs (%1.3fs CPU)\n", - i + 1, num_band0, loopTotalWallTime, loopTotalCPUTime); -/* #else */ -/* printf("*"); */ -/* if (i == (num_band0 - 1)) { */ -/* printf("\n"); */ -/* } */ -/* if ((i % 20) == 0 && i != 0) { */ -/* printf("\n"); */ -/* } */ + printf(" %1.3fs (%1.3fs CPU)\n", loopTotalWallTime, loopTotalCPUTime); #endif - } else { - for (j = 0; j < num_band * num_band; j++) { - fc3_normal_squared[i * num_band * num_band + j] = 0; - } - -#ifdef MEASURE_R2N - printf(" Band index %d/%d skipped due to frequency cutoff...\n", - i + 1, num_band0); -#endif - - } - } + free(n); } static double get_fc3_sum diff --git a/c/anharmonic_h/phonon3_h/real_to_reciprocal.h b/c/anharmonic_h/phonon3_h/real_to_reciprocal.h index de2c84b7..ad482f11 100644 --- a/c/anharmonic_h/phonon3_h/real_to_reciprocal.h +++ b/c/anharmonic_h/phonon3_h/real_to_reciprocal.h @@ -44,6 +44,7 @@ void real_to_reciprocal(lapack_complex_double *fc3_reciprocal, const Darray *shortest_vectors, const Iarray *multiplicity, const int *p2s_map, - const int *s2p_map); + const int *s2p_map, + const int openmp_at_bands); #endif diff --git a/scripts/phono3py b/scripts/phono3py index d802b8d9..d9c32bcb 100755 --- a/scripts/phono3py +++ b/scripts/phono3py @@ -644,6 +644,7 @@ elif run_mode == "conductivity-RTA" or run_mode == "conductivity-LBTE": is_kappa_star=settings.get_is_kappa_star(), gv_delta_q=settings.get_group_velocity_delta_q(), run_with_g=settings.get_run_with_g(), + with_g_zero=settings.get_with_g_zero(), pinv_cutoff=settings.get_pinv_cutoff(), write_gamma=settings.get_write_gamma(), read_gamma=settings.get_read_gamma(),