From c763d117c8262f9d7aa2822dd1cc511e34065e9d Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sun, 24 Dec 2023 19:09:33 +0900 Subject: [PATCH 1/2] fc3 transform with r0-average and all-shortest atom triplets --- c/_phono3py.c | 28 +++-- c/phono3py.c | 13 ++- c/phono3py.h | 11 +- c/real_to_reciprocal.c | 48 ++++++-- c/real_to_reciprocal.h | 1 + phono3py/conductivity/direct_solution.py | 4 +- phono3py/conductivity/rta.py | 14 +-- phono3py/phonon3/interaction.py | 143 +++++++++++++++-------- 8 files changed, 174 insertions(+), 88 deletions(-) diff --git a/c/_phono3py.c b/c/_phono3py.c index 33f396fc..e0d5d445 100644 --- a/c/_phono3py.c +++ b/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,6 +294,7 @@ 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; @@ -300,12 +302,12 @@ static PyObject *py_get_interaction(PyObject *self, PyObject *args) { 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; diff --git a/c/phono3py.c b/c/phono3py.c index 06deca27..133fd646 100644 --- a/c/phono3py.c +++ b/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, diff --git a/c/phono3py.h b/c/phono3py.h index 8c611543..febbd367 100644 --- a/c/phono3py.h +++ b/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], diff --git a/c/real_to_reciprocal.c b/c/real_to_reciprocal.c index 73fc4a6a..90a75f33 100644 --- a/c/real_to_reciprocal.c +++ b/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,6 +204,12 @@ 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, @@ -207,11 +217,27 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem, 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; diff --git a/c/real_to_reciprocal.h b/c/real_to_reciprocal.h index 721cbbdf..d8bef345 100644 --- a/c/real_to_reciprocal.h +++ b/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, diff --git a/phono3py/conductivity/direct_solution.py b/phono3py/conductivity/direct_solution.py index 02544964..132a4fd6 100644 --- a/phono3py/conductivity/direct_solution.py +++ b/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 diff --git a/phono3py/conductivity/rta.py b/phono3py/conductivity/rta.py index 204d2231..be3fd1f0 100644 --- a/phono3py/conductivity/rta.py +++ b/phono3py/conductivity/rta.py @@ -115,7 +115,7 @@ def __init__( log_level=log_level, ) - self._use_const_ave_pp = self._pp.get_constant_averaged_interaction() + 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 @@ -262,10 +262,10 @@ def _set_gamma_at_sigmas(self, i): if self._log_level: print( "Constant ph-ph interaction of %6.3e is used." - % self._pp.get_constant_averaged_interaction() + % self._pp.constant_averaged_interaction ) self._collision.run_interaction() - self._averaged_pp_interaction[i] = self._pp.get_averaged_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.") @@ -274,9 +274,7 @@ def _set_gamma_at_sigmas(self, i): 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.get_averaged_interaction() + self._averaged_pp_interaction[i] = self._pp.averaged_interaction # Number of triplets depends on q-point. # So this is allocated each time. @@ -386,6 +384,7 @@ class instance. 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, ) @@ -416,11 +415,12 @@ class instance. 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.get_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] diff --git a/phono3py/phonon3/interaction.py b/phono3py/phonon3/interaction.py index 9579e99b..f0ec2aaa 100644 --- a/phono3py/phonon3/interaction.py +++ b/phono3py/phonon3/interaction.py @@ -32,13 +32,14 @@ # 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 __future__ import annotations import warnings from collections.abc import Sequence from typing import Literal, Optional, Union import numpy as np -from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix +from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, get_dynamical_matrix from phonopy.structure.cells import ( Primitive, compute_all_sg_permutations, @@ -211,7 +212,7 @@ def run( ) @property - def interaction_strength(self): + def interaction_strength(self) -> np.ndarray: """Return ph-ph interaction strength. Returns @@ -233,7 +234,7 @@ def get_interaction_strength(self): return self.interaction_strength @property - def mesh_numbers(self): + def mesh_numbers(self) -> np.ndarray: """Return mesh numbers. Returns @@ -254,12 +255,12 @@ def get_mesh_numbers(self): return self.mesh_numbers @property - def is_mesh_symmetry(self): + def is_mesh_symmetry(self) -> bool: """Whether symmetry of grid is utilized or not.""" return self._is_mesh_symmetry @property - def fc3(self): + def fc3(self) -> np.ndarray: """Return fc3.""" return self._fc3 @@ -272,7 +273,7 @@ def get_fc3(self): return self.fc3 @property - def dynamical_matrix(self): + def dynamical_matrix(self) -> Optional[DynamicalMatrix]: """Return DynamicalMatrix class instance.""" return self._dm @@ -286,7 +287,7 @@ def get_dynamical_matrix(self): return self.dynamical_matrix @property - def primitive(self): + def primitive(self) -> Primitive: """Return Primitive class instance.""" return self._primitive @@ -300,11 +301,13 @@ def get_primitive(self): return self.primitive @property - def primitive_symmetry(self): + def primitive_symmetry(self) -> Symmetry: """Return Symmetry class instance of primitive cell.""" return self._primitive_symmetry - def get_triplets_at_q(self): + def get_triplets_at_q( + self, + ) -> tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray): """Return grid point triplets information. triplets_at_q is in BZ-grid. @@ -321,12 +324,12 @@ def get_triplets_at_q(self): ) @property - def bz_grid(self): + def bz_grid(self) -> BZGrid: """Return BZGrid class instance.""" return self._bz_grid @property - def band_indices(self): + def band_indices(self) -> np.ndarray: """Return band indices. Returns @@ -347,12 +350,12 @@ def get_band_indices(self): return self.band_indices @property - def nac_params(self): + def nac_params(self) -> dict: """Return NAC params.""" return self._nac_params @property - def nac_q_direction(self): + def nac_q_direction(self) -> Optional[np.ndarray]: """Return q-direction used for NAC at q->0. Direction of q-vector watching from Gamma point used for @@ -391,7 +394,7 @@ def set_nac_q_direction(self, nac_q_direction=None): self.nac_q_direction = nac_q_direction @property - def zero_value_positions(self): + def zero_value_positions(self) -> Optional[np.ndarray]: """Return zero ph-ph interaction elements information. Returns @@ -410,7 +413,7 @@ def get_zero_value_positions(self): ) return self.zero_value_positions - def get_phonons(self): + def get_phonons(self) -> tuple(np.ndarray, np.ndarray, np.ndarray): """Return phonons on grid. Returns @@ -431,7 +434,7 @@ def get_phonons(self): return self._frequencies, self._eigenvectors, self._phonon_done @property - def frequency_factor_to_THz(self): + def frequency_factor_to_THz(self) -> float: """Return phonon frequency conversion factor to THz.""" return self._frequency_factor_to_THz @@ -445,7 +448,7 @@ def get_frequency_factor_to_THz(self): return self.frequency_factor_to_THz @property - def lapack_zheev_uplo(self): + def lapack_zheev_uplo(self) -> Literal["L", "U"]: """Return U or L for lapack zheev solver.""" return self._lapack_zheev_uplo @@ -459,7 +462,7 @@ def get_lapack_zheev_uplo(self): return self.lapack_zheev_uplo @property - def cutoff_frequency(self): + def cutoff_frequency(self) -> float: """Return cutoff phonon frequency to judge imaginary phonon.""" return self._cutoff_frequency @@ -473,17 +476,17 @@ def get_cutoff_frequency(self): return self.cutoff_frequency @property - def openmp_per_triplets(self): + def openmp_per_triplets(self) -> bool: """Return whether OpenMP distribution over triplets or bands.""" return self._openmp_per_triplets @property - def symmetrize_fc3q(self): + def symmetrize_fc3q(self) -> bool: """Return boolean of symmetrize_fc3q.""" return self._symmetrize_fc3q @property - def make_r0_average(self): + def make_r0_average(self) -> bool: """Return boolean of make_r0_average. This flag is used to activate averaging of fc3 transformation @@ -494,7 +497,20 @@ def make_r0_average(self): """ return self._make_r0_average - def get_averaged_interaction(self): + @property + def all_shortest(self) -> np.ndarray: + """Return boolean of make_r0_average. + + This flag is used to activate averaging of fc3 transformation + from real space to reciprocal space around three atoms. With False, + it is done at the first atom. With True, it is done at three atoms + and averaged. + + """ + return self._all_shortest + + @property + def averaged_interaction(self) -> np.ndarray: """Return sum over phonon triplets of interaction strength. See Eq.(21) of PRB 91, 094306 (2015) @@ -506,18 +522,49 @@ def get_averaged_interaction(self): v_sum = np.dot(w, v.sum(axis=2).sum(axis=2)) return v_sum / np.prod(v.shape[2:]) - def get_primitive_and_supercell_correspondence(self): + def get_averaged_interaction(self): + """Return sum over phonon triplets of interaction strength.""" + warnings.warn( + "Use attribute, Interaction.averaged_interaction " + "instead of Interaction.get_averaged_interaction().", + DeprecationWarning, + ) + return self.averaged_interaction + + def get_primitive_and_supercell_correspondence( + self, + ) -> tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray): """Return atomic pair information.""" return (self._svecs, self._multi, self._p2s, self._s2p, self._masses) - def get_unit_conversion_factor(self): + @property + def unit_conversion_factor(self) -> float: """Return unit conversion factor.""" return self._unit_conversion - def get_constant_averaged_interaction(self): + def get_unit_conversion_factor(self): + """Return unit conversion factor.""" + warnings.warn( + "Use attribute, Interaction.unit_conversion_factor " + "instead of Interaction.get_unit_conversion_factor().", + DeprecationWarning, + ) + return self.unit_conversion_factor + + @property + def constant_averaged_interaction(self) -> float: """Return constant averaged interaction.""" return self._constant_averaged_interaction + def get_constant_averaged_interaction(self): + """Return constant averaged interaction.""" + warnings.warn( + "Use attribute, Interaction.constant_averaged_interaction " + "instead of Interaction.get_constant_averaged_interaction().", + DeprecationWarning, + ) + return self.constant_averaged_interaction + def set_interaction_strength(self, pp_strength, g_zero=None): """Set interaction strength.""" self._interaction_strength = pp_strength @@ -915,6 +962,7 @@ def _run_c(self, g_zero): self._band_indices, self._symmetrize_fc3q * 1, self._make_r0_average * 1, + self._all_shortest, self._cutoff_frequency, openmp_per_triplets * 1, ) @@ -1007,30 +1055,27 @@ def _get_all_shortest(self): s2pp_map = [self._primitive.p2p_map[i] for i in self._s2p] lattice = self._primitive.cell - for i_patom in range(n_patom): - for j_atom in range(n_satom): - j_patom = s2pp_map[j_atom] - i_perm = np.where(perms[:, j_atom] == self._p2s[j_patom])[0][0] - for k_atom in range(n_satom): - initial_vec = ( - svecs[multi[k_atom, i_patom, 1]] - - svecs[multi[j_atom, i_patom, 1]] - ) - d_jk_shortest = np.linalg.norm(initial_vec @ lattice) - for j_m, k_m in np.ndindex( - (multi[j_atom, i_patom, 0], multi[k_atom, i_patom, 0]) - ): - vec_ij = svecs[multi[j_atom, i_patom, 1] + j_m] - vec_ik = svecs[multi[k_atom, i_patom, 1] + k_m] - d_jk_attempt = np.linalg.norm((vec_ik - vec_ij) @ lattice) - if d_jk_attempt < d_jk_shortest: - d_jk_shortest = d_jk_attempt - k_atom_mapped = perms[i_perm, k_atom] - d_jk_mapped = np.linalg.norm( - svecs[multi[k_atom_mapped, j_patom, 1]] @ lattice - ) - if abs(d_jk_mapped - d_jk_shortest) < self._symprec: - self._all_shortest[i_patom, j_atom, k_atom] = 1 + for i_patom, j_atom in np.ndindex((n_patom, n_satom)): + if multi[j_atom, i_patom, 0] > 1: + continue + j_patom = s2pp_map[j_atom] + i_perm = np.where(perms[:, j_atom] == self._p2s[j_patom])[0] + assert len(i_perm) == 1 + for k_atom in range(n_satom): + if multi[k_atom, i_patom, 0] > 1: + continue + k_atom_mapped = perms[i_perm[0], k_atom] + if multi[k_atom_mapped, j_patom, 0] > 1: + continue + vec_jk = ( + svecs[multi[k_atom, i_patom, 1]] - svecs[multi[j_atom, i_patom, 1]] + ) + d_jk = np.linalg.norm(vec_jk @ lattice) + d_jk_mapped = np.linalg.norm( + svecs[multi[k_atom_mapped, j_patom, 1]] @ lattice + ) + if abs(d_jk_mapped - d_jk) < self._symprec: + self._all_shortest[i_patom, j_atom, k_atom] = 1 def all_bands_exist(interaction: Interaction): From fb6642ab94c4791a38a3e9cd1b5d48dc9bc9424c Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sun, 24 Dec 2023 19:30:05 +0900 Subject: [PATCH 2/2] Fix test_get_all_shortest --- test/phonon3/test_interaction.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/phonon3/test_interaction.py b/test/phonon3/test_interaction.py index 22f61b26..3ea907ac 100644 --- a/test/phonon3/test_interaction.py +++ b/test/phonon3/test_interaction.py @@ -314,16 +314,15 @@ def test_get_all_shortest(aln_lda: Phono3py): svecs, multi, _, _, _ = itr.get_primitive_and_supercell_correspondence() n_satom, n_patom, _ = multi.shape for i, j, k in np.ndindex((n_patom, n_satom, n_satom)): - d_jk_shortest = np.linalg.norm(s_svecs[s_multi[j, k, 1]] @ s_lattice) is_found = 0 - for m_j, m_k in np.ndindex((multi[j, i, 0], multi[k, i, 0])): - vec_ij = svecs[multi[j, i, 1] + m_j] - vec_ik = svecs[multi[k, i, 1] + m_k] + if multi[j, i, 0] == 1 and multi[k, i, 0] == 1 and s_multi[j, k, 0] == 1: + d_jk_shortest = np.linalg.norm(s_svecs[s_multi[j, k, 1]] @ s_lattice) + vec_ij = svecs[multi[j, i, 1]] + vec_ik = svecs[multi[k, i, 1]] vec_jk = vec_ik - vec_ij d_jk = np.linalg.norm(vec_jk @ p_lattice) if abs(d_jk - d_jk_shortest) < ph3.symmetry.tolerance: is_found = 1 - break assert shortests[i, j, k] == is_found