From 4593caa9d2dfe3f1ee4d1054cb18d9ae46e6319f Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Tue, 19 Mar 2024 15:00:04 +0900 Subject: [PATCH 01/16] Follow recent change of phonopy --- phono3py/other/isotope.py | 1 - phono3py/phonon3/gruneisen.py | 1 - phono3py/phonon3/interaction.py | 1 - phono3py/phonon3/joint_dos.py | 1 - test/phonon/test_velocity_operator.py | 2 +- 5 files changed, 1 insertion(+), 5 deletions(-) diff --git a/phono3py/other/isotope.py b/phono3py/other/isotope.py index a786c08f..f929e2db 100644 --- a/phono3py/other/isotope.py +++ b/phono3py/other/isotope.py @@ -256,7 +256,6 @@ def init_dynamical_matrix( nac_params=nac_params, frequency_scale_factor=frequency_scale_factor, decimals=decimals, - symprec=self._symprec, ) def set_nac_q_direction(self, nac_q_direction=None): diff --git a/phono3py/phonon3/gruneisen.py b/phono3py/phonon3/gruneisen.py index c803d22b..bc47a3f9 100644 --- a/phono3py/phonon3/gruneisen.py +++ b/phono3py/phonon3/gruneisen.py @@ -149,7 +149,6 @@ def __init__( self._scell, self._pcell, nac_params=nac_params, - symprec=self._symprec, ) self._nac_q_direction = nac_q_direction diff --git a/phono3py/phonon3/interaction.py b/phono3py/phonon3/interaction.py index 581250d9..575088e6 100644 --- a/phono3py/phonon3/interaction.py +++ b/phono3py/phonon3/interaction.py @@ -671,7 +671,6 @@ def init_dynamical_matrix( nac_params=nac_params, frequency_scale_factor=self._frequency_scale_factor, decimals=decimals, - symprec=self._symprec, ) self._allocate_phonon() diff --git a/phono3py/phonon3/joint_dos.py b/phono3py/phonon3/joint_dos.py index 69fe9007..66d2e4da 100644 --- a/phono3py/phonon3/joint_dos.py +++ b/phono3py/phonon3/joint_dos.py @@ -417,7 +417,6 @@ def _init_dynamical_matrix(self): self._primitive, nac_params=self._nac_params, frequency_scale_factor=self._frequency_scale_factor, - symprec=self._symprec, ) self._allocate_phonons() diff --git a/test/phonon/test_velocity_operator.py b/test/phonon/test_velocity_operator.py index 1a51e700..320f2dd7 100644 --- a/test/phonon/test_velocity_operator.py +++ b/test/phonon/test_velocity_operator.py @@ -189,7 +189,7 @@ def test_gv_operator_nacl(ph_nacl: Phonopy): ) # we chose an 'ugly' q-point because we want to avoid degeneracies. - ph_nacl.dynamical_matrix.run([[0.1, 0.22, 0.33]]) + ph_nacl.dynamical_matrix.run([0.1, 0.22, 0.33]) dm = ph_nacl.dynamical_matrix.dynamical_matrix eigvals, eigvecs = np.linalg.eigh(dm) From f18f952a4814194e1c08ba192b8fc71febfc43b0 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Wed, 20 Mar 2024 09:41:29 +0900 Subject: [PATCH 02/16] Set version 2.10.0 --- doc/changelog.md | 4 ++++ doc/conf.py | 4 ++-- doc/install.md | 2 +- phono3py/version.py | 2 +- requirements.txt | 2 +- setup.py | 2 +- 6 files changed, 10 insertions(+), 6 deletions(-) diff --git a/doc/changelog.md b/doc/changelog.md index a2cfd955..c6911716 100644 --- a/doc/changelog.md +++ b/doc/changelog.md @@ -2,6 +2,10 @@ # Change Log +## Mar-20-2024: Version 2.10.0 + +- Maintenance release + ## Feb-2-2024: Version 2.9.2 - `boundary_mfp` value is stored in `kappa-*.hdf5` file when it is specified. diff --git a/doc/conf.py b/doc/conf.py index cae6e3a1..140fdb53 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -58,9 +58,9 @@ # built documents. # # The short X.Y version. -version = "2.9" +version = "2.10" # The full version, including alpha/beta/rc tags. -release = "2.9.2" +release = "2.10.0" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/install.md b/doc/install.md index 834c6314..16d64e77 100644 --- a/doc/install.md +++ b/doc/install.md @@ -161,7 +161,7 @@ wrong python libraries can be imported. For macOS ARM64 system, currently only openblas can be chosen: ```bash - % conda install numpy scipy h5py pyyaml matplotlib-base c-compiler spglib cmake openblas="0.3.18" + % conda install numpy scipy h5py pyyaml matplotlib-base c-compiler spglib cmake openblas ``` Note that using hdf5 files on NFS mounted file system, you may have to disable diff --git a/phono3py/version.py b/phono3py/version.py index 434dd992..50d3a786 100644 --- a/phono3py/version.py +++ b/phono3py/version.py @@ -34,4 +34,4 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -__version__ = "2.9.2" +__version__ = "2.10.0" diff --git a/requirements.txt b/requirements.txt index 4e13b8dc..ed682251 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,4 @@ numpy >= 1.17.0 PyYAML >= 5.3 matplotlib >= 2.2.2 h5py >= 3.0 -phonopy >=2.21,<2.22 +phonopy >=2.22,<2.23 diff --git a/setup.py b/setup.py index 7f0c1780..d221f990 100644 --- a/setup.py +++ b/setup.py @@ -320,7 +320,7 @@ def main(build_dir): "matplotlib>=2.2.2", "h5py>=3.0", "spglib>=2.0", - "phonopy>=2.21,<2.22", + "phonopy>=2.22,<2.23", ], provides=["phono3py"], scripts=scripts_phono3py, From a5a5615698d849c9e761e0269afe4db11c9df997 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 25 Mar 2024 22:32:21 +0000 Subject: [PATCH 03/16] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/psf/black: 24.2.0 → 24.3.0](https://github.com/psf/black/compare/24.2.0...24.3.0) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 59c4c9a3..d0c337a6 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -19,7 +19,7 @@ repos: - "--ignore=E203,W503" - repo: https://github.com/psf/black - rev: 24.2.0 + rev: 24.3.0 hooks: - id: black args: From b6ffd73720b4ff6ca4b19bdb0e9b9b9e0ff86ae0 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 8 Apr 2024 23:20:58 +0000 Subject: [PATCH 04/16] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/pre-commit/pre-commit-hooks: v4.5.0 → v4.6.0](https://github.com/pre-commit/pre-commit-hooks/compare/v4.5.0...v4.6.0) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index d0c337a6..4bc66156 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -2,7 +2,7 @@ # See https://pre-commit.com/hooks.html for more hooks repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.5.0 + rev: v4.6.0 hooks: - id: trailing-whitespace - id: end-of-file-fixer From 15e76c5f82066a2843e9643855061addc2fc8122 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 15 Apr 2024 23:05:54 +0000 Subject: [PATCH 05/16] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/psf/black: 24.3.0 → 24.4.0](https://github.com/psf/black/compare/24.3.0...24.4.0) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 4bc66156..ab24f699 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -19,7 +19,7 @@ repos: - "--ignore=E203,W503" - repo: https://github.com/psf/black - rev: 24.3.0 + rev: 24.4.0 hooks: - id: black args: From f28c6b17414bd1223c60aece0c16713c45520c42 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Wed, 17 Apr 2024 08:47:33 +0900 Subject: [PATCH 06/16] Minor refactoring of reciprocal_to_normal.c --- c/reciprocal_to_normal.c | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/c/reciprocal_to_normal.c b/c/reciprocal_to_normal.c index e55c6908..bf30e175 100644 --- a/c/reciprocal_to_normal.c +++ b/c/reciprocal_to_normal.c @@ -158,9 +158,9 @@ static double get_fc3_sum(const lapack_complex_double *e0, const lapack_complex_double *e2, const lapack_complex_double *fc3_reciprocal, const long num_band) { - long i, j, k; + long i, j, jk; double sum_real, sum_imag; - lapack_complex_double e_012, e_012_fc3, *e_12_cache; + lapack_complex_double e_012_fc3, fc3_i_e_12, *e_12_cache; const lapack_complex_double *fc3_i; e_12_cache = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) * @@ -177,19 +177,15 @@ static double get_fc3_sum(const lapack_complex_double *e0, for (i = 0; i < num_band; i++) { fc3_i = fc3_reciprocal + i * num_band * num_band; - for (j = 0; j < num_band; j++) { - for (k = 0; k < num_band; k++) { - e_012 = - phonoc_complex_prod(e0[i], e_12_cache[j * num_band + k]); - e_012_fc3 = phonoc_complex_prod(e_012, fc3_i[j * num_band + k]); - sum_real += lapack_complex_double_real(e_012_fc3); - sum_imag += lapack_complex_double_imag(e_012_fc3); - } + for (jk = 0; jk < num_band * num_band; jk++) { + fc3_i_e_12 = phonoc_complex_prod(fc3_i[jk], e_12_cache[jk]); + e_012_fc3 = phonoc_complex_prod(e0[i], fc3_i_e_12); + sum_real += lapack_complex_double_real(e_012_fc3); + sum_imag += lapack_complex_double_imag(e_012_fc3); } } free(e_12_cache); e_12_cache = NULL; - return (sum_real * sum_real + sum_imag * sum_imag); } From 16f3087267c9e6c8e0ff0c3028c2220b4ceabda1 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Wed, 17 Apr 2024 09:08:29 +0900 Subject: [PATCH 07/16] Memorize phase factors in inner loop for better performance in calculation --- c/real_to_reciprocal.c | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/c/real_to_reciprocal.c b/c/real_to_reciprocal.c index db54dd53..7adde156 100644 --- a/c/real_to_reciprocal.c +++ b/c/real_to_reciprocal.c @@ -233,22 +233,31 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem, 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; + lapack_complex_double phase_factor, phase_factor1, *phase_factor2; double fc3_rec_real[27], fc3_rec_imag[27]; + num_satom = atom_triplets->multi_dims[0]; + + phase_factor2 = (lapack_complex_double *)malloc( + sizeof(lapack_complex_double) * num_satom); + for (i = 0; i < 27; i++) { fc3_rec_real[i] = 0; fc3_rec_imag[i] = 0; } - num_satom = atom_triplets->multi_dims[0]; - if (is_compact_fc3) { i = pi0; } else { i = atom_triplets->p2s_map[pi0]; } + for (j = 0; j < num_satom; j++) { + adrs_vec2 = j * atom_triplets->multi_dims[1] + pi0; + phase_factor2[j] = get_phase_factor( + q2, atom_triplets->svecs, atom_triplets->multiplicity[adrs_vec2]); + } + for (j = 0; j < num_satom; j++) { if (atom_triplets->s2p_map[j] != atom_triplets->p2s_map[pi1]) { continue; @@ -267,13 +276,9 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem, 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); + phase_factor = phonoc_complex_prod(phase_factor1, phase_factor2[k]); if ((leg_index == 1) && (atom_triplets->all_shortest[pi0 * num_satom * num_satom + @@ -303,6 +308,9 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem, fc3_rec_elem[i] = lapack_make_complex_double(fc3_rec_real[i], fc3_rec_imag[i]); } + + free(phase_factor2); + phase_factor2 = NULL; } // This function doesn't need to think about position + From 293ede7719cf55c0482bff7fb54db971fab1f33a Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Wed, 17 Apr 2024 15:36:26 +0900 Subject: [PATCH 08/16] Examine BLAS for ph-ph interaction calculation --- c/reciprocal_to_normal.c | 44 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/c/reciprocal_to_normal.c b/c/reciprocal_to_normal.c index bf30e175..a17e7cf5 100644 --- a/c/reciprocal_to_normal.c +++ b/c/reciprocal_to_normal.c @@ -34,6 +34,7 @@ #include "reciprocal_to_normal.h" +#include #include #include @@ -49,7 +50,11 @@ static double get_fc3_sum(const lapack_complex_double *e0, const lapack_complex_double *e2, const lapack_complex_double *fc3_reciprocal, const long num_band); - +static double get_fc3_sum_blas(const lapack_complex_double *e0, + const lapack_complex_double *e1, + const lapack_complex_double *e2, + const lapack_complex_double *fc3_reciprocal, + const long num_band); void reciprocal_to_normal_squared( double *fc3_normal_squared, const long (*g_pos)[4], const long num_g_pos, const lapack_complex_double *fc3_reciprocal, const double *freqs0, @@ -189,3 +194,40 @@ static double get_fc3_sum(const lapack_complex_double *e0, e_12_cache = NULL; return (sum_real * sum_real + sum_imag * sum_imag); } + +static double get_fc3_sum_blas(const lapack_complex_double *e0, + const lapack_complex_double *e1, + const lapack_complex_double *e2, + const lapack_complex_double *fc3_reciprocal, + const long num_band) { + long i, j; + lapack_complex_double *fc3_e12, *e_12, zero, one, retval; + const lapack_complex_double *fc3_i; + + e_12 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) * + num_band * num_band); + fc3_e12 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) * + num_band); + zero = lapack_make_complex_double(0, 0); + one = lapack_make_complex_double(1, 0); + + for (i = 0; i < num_band; i++) { + cblas_zcopy(num_band, e2, 1, e_12 + i * num_band, 1); + cblas_zscal(num_band, e1 + i, e_12 + i * num_band, 1); + } + + cblas_zgemv(CblasRowMajor, CblasNoTrans, num_band, num_band * num_band, + &one, fc3_reciprocal, num_band * num_band, e_12, 1, &zero, + fc3_e12, 1); + cblas_zdotu_sub(num_band, e0, 1, fc3_e12, 1, &retval); + + free(e_12); + e_12 = NULL; + free(fc3_e12); + fc3_e12 = NULL; + + return lapack_complex_double_real(retval) * + lapack_complex_double_real(retval) + + lapack_complex_double_imag(retval) * + lapack_complex_double_imag(retval); +} From a46e234729d28b275e67a5fc3c3e4370655253b7 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Wed, 17 Apr 2024 16:05:11 +0900 Subject: [PATCH 09/16] cblas header for mkl --- c/reciprocal_to_normal.c | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/c/reciprocal_to_normal.c b/c/reciprocal_to_normal.c index a17e7cf5..c3f39bef 100644 --- a/c/reciprocal_to_normal.c +++ b/c/reciprocal_to_normal.c @@ -34,7 +34,11 @@ #include "reciprocal_to_normal.h" +#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) +#include +#else #include +#endif #include #include From 44a404abe43144e9964146a6a9d91019d5657786 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Thu, 18 Apr 2024 12:55:38 +0900 Subject: [PATCH 10/16] Replace openmp_(at|per)_bands by (not openmp_per_triplets) in C --- c/_phono3py.c | 2 +- c/gridsys.c | 10 +- c/gridsys.h | 2 +- c/interaction.c | 24 ++--- c/interaction.h | 2 +- c/phono3py.c | 10 +- c/phono3py.h | 2 +- c/pp_collision.c | 4 +- c/real_to_reciprocal.c | 202 +++++++++++++++++++++++---------------- c/real_to_reciprocal.h | 2 +- c/reciprocal_to_normal.c | 9 +- c/reciprocal_to_normal.h | 2 +- c/triplet.c | 4 +- c/triplet.h | 2 +- c/triplet_iw.c | 8 +- c/triplet_iw.h | 4 +- 16 files changed, 161 insertions(+), 128 deletions(-) diff --git a/c/_phono3py.c b/c/_phono3py.c index e0d5d445..27388cb5 100644 --- a/c/_phono3py.c +++ b/c/_phono3py.c @@ -1419,7 +1419,7 @@ static PyObject *py_get_triplets_integration_weights(PyObject *self, ph3py_get_integration_weight( iw, iw_zero, frequency_points, num_band0, relative_grid_address, D_diag, triplets, num_triplets, bz_grid_addresses, bz_map, bz_grid_type, - frequencies1, num_band1, frequencies2, num_band2, tp_type, 1, 0); + frequencies1, num_band1, frequencies2, num_band2, tp_type, 1); Py_RETURN_NONE; } diff --git a/c/gridsys.c b/c/gridsys.c index aa10eb68..4c8aa645 100644 --- a/c/gridsys.c +++ b/c/gridsys.c @@ -266,7 +266,7 @@ long gridsys_get_integration_weight( const long (*bz_grid_addresses)[3], const long *bz_map, const long bz_grid_type, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets, const long openmp_per_bands) { + const long openmp_per_triplets) { ConstBZGrid *bzgrid; long i; @@ -282,10 +282,10 @@ long gridsys_get_integration_weight( bzgrid->D_diag[i] = D_diag[i]; } - tpl_get_integration_weight( - iw, iw_zero, frequency_points, num_band0, relative_grid_address, - triplets, num_triplets, bzgrid, frequencies1, num_band1, frequencies2, - num_band2, tp_type, openmp_per_triplets, openmp_per_bands); + tpl_get_integration_weight(iw, iw_zero, frequency_points, num_band0, + relative_grid_address, triplets, num_triplets, + bzgrid, frequencies1, num_band1, frequencies2, + num_band2, tp_type, openmp_per_triplets); free(bzgrid); bzgrid = NULL; diff --git a/c/gridsys.h b/c/gridsys.h index ebee0bdd..a7aa8431 100644 --- a/c/gridsys.h +++ b/c/gridsys.h @@ -310,7 +310,7 @@ long gridsys_get_integration_weight( const long (*bz_grid_addresses)[3], const long *bz_map, const long bz_grid_type, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets, const long openmp_per_bands); + const long openmp_per_triplets); void gridsys_get_integration_weight_with_sigma( double *iw, char *iw_zero, const double sigma, const double sigma_cutoff, const double *frequency_points, const long num_band0, diff --git a/c/interaction.c b/c/interaction.c index bce45e9b..cc7ba929 100644 --- a/c/interaction.c +++ b/c/interaction.c @@ -58,7 +58,7 @@ static void real_to_normal(double *fc3_normal_squared, const long (*g_pos)[4], const double *masses, const long *band_indices, const long num_band, const double cutoff_frequency, const long triplet_index, const long num_triplets, - const long openmp_at_bands); + const long openmp_per_triplets); static void real_to_normal_sym_q( double *fc3_normal_squared, const long (*g_pos)[4], const long num_g_pos, double *const freqs[3], lapack_complex_double *const eigvecs[3], @@ -67,7 +67,7 @@ static void real_to_normal_sym_q( const AtomTriplets *atom_triplets, const double *masses, const long *band_indices, const long num_band0, const long num_band, const double cutoff_frequency, const long triplet_index, - const long num_triplets, const long openmp_at_bands); + const long num_triplets, const long openmp_per_triplets); /* fc3_normal_squared[num_triplets, num_band0, num_band, num_band] */ void itr_get_interaction( @@ -101,7 +101,7 @@ void itr_get_interaction( g_pos, num_g_pos, frequencies->data, eigenvectors, triplets[i], bzgrid, fc3, is_compact_fc3, atom_triplets, masses, band_indices, symmetrize_fc3_q, cutoff_frequency, i, num_triplets, - 1 - openmp_per_triplets); + openmp_per_triplets); free(g_pos); g_pos = NULL; @@ -118,7 +118,7 @@ void itr_get_interaction_at_triplet( const double cutoff_frequency, const long triplet_index, /* only for print */ const long num_triplets, /* only for print */ - const long openmp_at_bands) { + const long openmp_per_triplets) { long j, k; double *freqs[3]; lapack_complex_double *eigvecs[3]; @@ -149,7 +149,7 @@ void itr_get_interaction_at_triplet( fc3_normal_squared, g_pos, num_g_pos, freqs, eigvecs, fc3, is_compact_fc3, q_vecs, /* q0, q1, q2 */ atom_triplets, masses, band_indices, num_band0, num_band, - cutoff_frequency, triplet_index, num_triplets, openmp_at_bands); + cutoff_frequency, triplet_index, num_triplets, openmp_per_triplets); for (j = 0; j < 3; j++) { free(freqs[j]); freqs[j] = NULL; @@ -167,7 +167,7 @@ void itr_get_interaction_at_triplet( is_compact_fc3, q_vecs, /* q0, q1, q2 */ atom_triplets, masses, band_indices, num_band, cutoff_frequency, triplet_index, num_triplets, - openmp_at_bands); + openmp_per_triplets); } } @@ -183,7 +183,7 @@ static void real_to_normal(double *fc3_normal_squared, const long (*g_pos)[4], const double *masses, const long *band_indices, const long num_band, const double cutoff_frequency, const long triplet_index, const long num_triplets, - const long openmp_at_bands) { + const long openmp_per_triplets) { lapack_complex_double *fc3_reciprocal; lapack_complex_double comp_zero; long i; @@ -195,10 +195,10 @@ static void real_to_normal(double *fc3_normal_squared, const long (*g_pos)[4], fc3_reciprocal[i] = comp_zero; } r2r_real_to_reciprocal(fc3_reciprocal, q_vecs, fc3, is_compact_fc3, - atom_triplets, openmp_at_bands); + atom_triplets, openmp_per_triplets); #ifdef MEASURE_R2N - if (openmp_at_bands && num_triplets > 0) { + if ((!openmp_per_triplets) && num_triplets > 0) { printf("At triplet %d/%d (# of bands=%d):\n", triplet_index, num_triplets, num_band0); } @@ -206,7 +206,7 @@ static void real_to_normal(double *fc3_normal_squared, const long (*g_pos)[4], 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_at_bands); + cutoff_frequency, openmp_per_triplets); free(fc3_reciprocal); fc3_reciprocal = NULL; @@ -220,7 +220,7 @@ static void real_to_normal_sym_q( const AtomTriplets *atom_triplets, const double *masses, const long *band_indices, const long num_band0, const long num_band, const double cutoff_frequency, const long triplet_index, - const long num_triplets, const long openmp_at_bands) { + const long num_triplets, const long openmp_per_triplets) { long i, j, k, l; long band_ex[3]; double q_vecs_ex[3][3]; @@ -246,7 +246,7 @@ static void real_to_normal_sym_q( 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_at_bands); + 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++) { diff --git a/c/interaction.h b/c/interaction.h index 2be35927..c345f742 100644 --- a/c/interaction.h +++ b/c/interaction.h @@ -57,6 +57,6 @@ void itr_get_interaction_at_triplet( const double cutoff_frequency, const long triplet_index, /* only for print */ const long num_triplets, /* only for print */ - const long openmp_at_bands); + const long openmp_per_triplets); #endif diff --git a/c/phono3py.c b/c/phono3py.c index 133fd646..e5cee178 100644 --- a/c/phono3py.c +++ b/c/phono3py.c @@ -435,7 +435,7 @@ long ph3py_get_integration_weight( const long (*bz_grid_addresses)[3], const long *bz_map, const long bz_grid_type, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets, const long openmp_per_bands) { + const long openmp_per_triplets) { ConstBZGrid *bzgrid; long i; @@ -451,10 +451,10 @@ long ph3py_get_integration_weight( bzgrid->D_diag[i] = D_diag[i]; } - tpl_get_integration_weight( - iw, iw_zero, frequency_points, num_band0, relative_grid_address, - triplets, num_triplets, bzgrid, frequencies1, num_band1, frequencies2, - num_band2, tp_type, openmp_per_triplets, openmp_per_bands); + tpl_get_integration_weight(iw, iw_zero, frequency_points, num_band0, + relative_grid_address, triplets, num_triplets, + bzgrid, frequencies1, num_band1, frequencies2, + num_band2, tp_type, openmp_per_triplets); free(bzgrid); bzgrid = NULL; diff --git a/c/phono3py.h b/c/phono3py.h index febbd367..1d9ad000 100644 --- a/c/phono3py.h +++ b/c/phono3py.h @@ -166,7 +166,7 @@ long ph3py_get_integration_weight( const long (*bz_grid_addresses)[3], const long *bz_map, const long bz_grid_type, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets, const long openmp_per_bands); + const long openmp_per_triplets); void ph3py_get_integration_weight_with_sigma( double *iw, char *iw_zero, const double sigma, const double sigma_cutoff, const double *frequency_points, const long num_band0, diff --git a/c/pp_collision.c b/c/pp_collision.c index 955ea38d..6258f083 100644 --- a/c/pp_collision.c +++ b/c/pp_collision.c @@ -109,7 +109,7 @@ void ppc_get_pp_collision( triplets[i], 1, bzgrid, frequencies, /* used as f1 */ num_band, frequencies, /* used as f2 */ - num_band, 2, 1 - openmp_per_triplets); + 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, @@ -179,7 +179,7 @@ void ppc_get_pp_collision_with_sigma( g_zero = (char *)malloc(sizeof(char) * num_band_prod); tpi_get_integration_weight_with_sigma( g, g_zero, sigma, cutoff, freqs_at_gp, num_band0, triplets[i], - const_adrs_shift, frequencies, num_band, 2, 0); + 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, diff --git a/c/real_to_reciprocal.c b/c/real_to_reciprocal.c index 7adde156..00f88068 100644 --- a/c/real_to_reciprocal.c +++ b/c/real_to_reciprocal.c @@ -43,25 +43,27 @@ #include "phonoc_array.h" #include "phonoc_const.h" -static void real_to_reciprocal_legacy(lapack_complex_double *fc3_reciprocal, - const double q_vecs[3][3], - const double *fc3, - const long is_compact_fc3, - const AtomTriplets *atom_triplets, - const long openmp_at_bands); -static void real_to_reciprocal_r0_average(lapack_complex_double *fc3_reciprocal, - const double q_vecs[3][3], - const double *fc3, - const long is_compact_fc3, - const AtomTriplets *atom_triplets, - const long openmp_at_bands); -static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem, - const double q1[3], const double q2[3], - const double *fc3, - const long is_compact_fc3, - const AtomTriplets *atom_triplets, - const long pi0, const long pi1, - const long pi2, const long leg_index); +static void real_to_reciprocal_legacy( + lapack_complex_double *fc3_reciprocal, + const lapack_complex_double *pre_phase_factors, + const lapack_complex_double *phase_factor1, + const lapack_complex_double *phase_factor2, const double *fc3, + const long is_compact_fc3, const AtomTriplets *atom_triplets, + const long openmp_per_triplets); +static void real_to_reciprocal_r0_average( + lapack_complex_double *fc3_reciprocal, + const lapack_complex_double *pre_phase_factors, + const lapack_complex_double *phase_factor0, + const lapack_complex_double *phase_factor1, + const lapack_complex_double *phase_factor2, const double *fc3, + const long is_compact_fc3, const AtomTriplets *atom_triplets, + const long openmp_per_triplets); +static void real_to_reciprocal_elements( + lapack_complex_double *fc3_rec_elem, + const lapack_complex_double *phase_factor1, + const lapack_complex_double *phase_factor2, const double *fc3, + const long is_compact_fc3, const AtomTriplets *atom_triplets, + const long pi0, const long pi1, 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]); @@ -76,13 +78,46 @@ void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal, const double q_vecs[3][3], const double *fc3, const long is_compact_fc3, const AtomTriplets *atom_triplets, - const long openmp_at_bands) { - long i, num_band; + const long openmp_per_triplets) { + long i, j, num_band, num_patom, num_satom, adrs_vec; + lapack_complex_double *pre_phase_factors, *phase_factor0, *phase_factor1, + *phase_factor2; + + num_patom = atom_triplets->multi_dims[1]; + num_satom = atom_triplets->multi_dims[0]; + + pre_phase_factors = (lapack_complex_double *)malloc( + sizeof(lapack_complex_double) * num_patom); + for (i = 0; i < num_patom; i++) { + pre_phase_factors[i] = get_pre_phase_factor(i, q_vecs, atom_triplets); + } + + phase_factor0 = (lapack_complex_double *)malloc( + sizeof(lapack_complex_double) * num_patom * num_satom); + phase_factor1 = (lapack_complex_double *)malloc( + sizeof(lapack_complex_double) * num_patom * num_satom); + phase_factor2 = (lapack_complex_double *)malloc( + sizeof(lapack_complex_double) * num_patom * num_satom); + for (i = 0; i < num_patom; i++) { + for (j = 0; j < num_satom; j++) { + adrs_vec = j * atom_triplets->multi_dims[1] + i; + phase_factor0[i * num_satom + j] = + get_phase_factor(q_vecs[0], atom_triplets->svecs, + atom_triplets->multiplicity[adrs_vec]); + phase_factor1[i * num_satom + j] = + get_phase_factor(q_vecs[1], atom_triplets->svecs, + atom_triplets->multiplicity[adrs_vec]); + phase_factor2[i * num_satom + j] = + get_phase_factor(q_vecs[2], atom_triplets->svecs, + atom_triplets->multiplicity[adrs_vec]); + } + } if (atom_triplets->make_r0_average) { - real_to_reciprocal_r0_average(fc3_reciprocal, q_vecs, fc3, - is_compact_fc3, atom_triplets, - openmp_at_bands); + real_to_reciprocal_r0_average(fc3_reciprocal, pre_phase_factors, + phase_factor0, phase_factor1, + phase_factor2, fc3, is_compact_fc3, + atom_triplets, openmp_per_triplets); num_band = atom_triplets->multi_dims[1] * 3; for (i = 0; i < num_band * num_band * num_band; i++) { fc3_reciprocal[i] = lapack_make_complex_double( @@ -90,42 +125,53 @@ void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal, lapack_complex_double_imag(fc3_reciprocal[i]) / 3); } } else { - real_to_reciprocal_legacy(fc3_reciprocal, q_vecs, fc3, is_compact_fc3, - atom_triplets, openmp_at_bands); + real_to_reciprocal_legacy( + fc3_reciprocal, pre_phase_factors, phase_factor1, phase_factor2, + fc3, is_compact_fc3, atom_triplets, openmp_per_triplets); } + + free(pre_phase_factors); + pre_phase_factors = NULL; + free(phase_factor0); + phase_factor1 = NULL; + free(phase_factor1); + phase_factor1 = NULL; + free(phase_factor2); + phase_factor2 = NULL; } -static void real_to_reciprocal_legacy(lapack_complex_double *fc3_reciprocal, - const double q_vecs[3][3], - const double *fc3, - const long is_compact_fc3, - const AtomTriplets *atom_triplets, - const long openmp_at_bands) { +static void real_to_reciprocal_legacy( + lapack_complex_double *fc3_reciprocal, + const lapack_complex_double *pre_phase_factors, + const lapack_complex_double *phase_factor1, + const lapack_complex_double *phase_factor2, const double *fc3, + const long is_compact_fc3, const AtomTriplets *atom_triplets, + const long openmp_per_triplets) { long i, j, k, l, m, n, ijk; - long num_patom, num_band; - lapack_complex_double pre_phase_factor, fc3_rec_elem[27], fc3_rec; + long num_patom, num_satom, num_band; + lapack_complex_double fc3_rec_elem[27], fc3_rec; num_patom = atom_triplets->multi_dims[1]; + num_satom = atom_triplets->multi_dims[0]; num_band = num_patom * 3; #ifdef _OPENMP -#pragma omp parallel for private(i, j, k, l, m, n, fc3_rec_elem, fc3_rec, \ - pre_phase_factor) if (openmp_at_bands) +#pragma omp parallel for private(i, j, k, l, m, n, fc3_rec_elem, \ + fc3_rec) if (!openmp_per_triplets) #endif for (ijk = 0; ijk < num_patom * num_patom * num_patom; ijk++) { 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, + real_to_reciprocal_elements(fc3_rec_elem, phase_factor1 + i * num_satom, + phase_factor2 + i * num_satom, fc3, is_compact_fc3, atom_triplets, i, j, k, 0); for (l = 0; l < 3; l++) { for (m = 0; m < 3; m++) { for (n = 0; n < 3; n++) { fc3_rec = phonoc_complex_prod( - fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factor); + fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factors[i]); fc3_reciprocal[(i * 3 + l) * num_band * num_band + (j * 3 + m) * num_band + k * 3 + n] = sum_lapack_complex_double( @@ -144,36 +190,39 @@ static void real_to_reciprocal_legacy(lapack_complex_double *fc3_reciprocal, // q-points in triplets used for summation implemented in // real_to_reciprocal_elements(). // --sym-fc3q makes them almost equivalent. -static void real_to_reciprocal_r0_average(lapack_complex_double *fc3_reciprocal, - const double q_vecs[3][3], - const double *fc3, - const long is_compact_fc3, - const AtomTriplets *atom_triplets, - const long openmp_at_bands) { +static void real_to_reciprocal_r0_average( + lapack_complex_double *fc3_reciprocal, + const lapack_complex_double *pre_phase_factors, + const lapack_complex_double *phase_factor0, + const lapack_complex_double *phase_factor1, + const lapack_complex_double *phase_factor2, const double *fc3, + const long is_compact_fc3, const AtomTriplets *atom_triplets, + const long openmp_per_triplets) { long i, j, k, l, m, n, ijk; - long num_patom, num_band; - lapack_complex_double pre_phase_factor, fc3_rec_elem[27], fc3_rec; + long num_patom, num_satom, num_band; + lapack_complex_double fc3_rec_elem[27], fc3_rec; num_patom = atom_triplets->multi_dims[1]; + num_satom = atom_triplets->multi_dims[0]; num_band = num_patom * 3; #ifdef _OPENMP -#pragma omp parallel for private(i, j, k, l, m, n, fc3_rec_elem, fc3_rec, \ - pre_phase_factor) if (openmp_at_bands) +#pragma omp parallel for private(i, j, k, l, m, n, fc3_rec_elem, \ + fc3_rec) if (!openmp_per_triplets) #endif for (ijk = 0; ijk < num_patom * num_patom * num_patom; ijk++) { 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, + real_to_reciprocal_elements(fc3_rec_elem, phase_factor1 + i * num_satom, + phase_factor2 + i * num_satom, fc3, is_compact_fc3, atom_triplets, i, j, k, 1); for (l = 0; l < 3; l++) { for (m = 0; m < 3; m++) { for (n = 0; n < 3; n++) { fc3_rec = phonoc_complex_prod( - fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factor); + fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factors[i]); fc3_reciprocal[(i * 3 + l) * num_band * num_band + (j * 3 + m) * num_band + k * 3 + n] = sum_lapack_complex_double( @@ -185,14 +234,14 @@ static void real_to_reciprocal_r0_average(lapack_complex_double *fc3_reciprocal, } // fc3_rec is stored in a way swapping jm <-> il. - pre_phase_factor = get_pre_phase_factor(j, q_vecs, atom_triplets); - real_to_reciprocal_elements(fc3_rec_elem, q_vecs[0], q_vecs[2], fc3, + real_to_reciprocal_elements(fc3_rec_elem, phase_factor0 + j * num_satom, + phase_factor2 + j * num_satom, fc3, is_compact_fc3, atom_triplets, j, i, k, 2); for (l = 0; l < 3; l++) { for (m = 0; m < 3; m++) { for (n = 0; n < 3; n++) { fc3_rec = phonoc_complex_prod( - fc3_rec_elem[m * 9 + l * 3 + n], pre_phase_factor); + fc3_rec_elem[m * 9 + l * 3 + n], pre_phase_factors[j]); fc3_reciprocal[(i * 3 + l) * num_band * num_band + (j * 3 + m) * num_band + k * 3 + n] = sum_lapack_complex_double( @@ -204,14 +253,14 @@ static void real_to_reciprocal_r0_average(lapack_complex_double *fc3_reciprocal, } // fc3_rec is stored in a way swapping kn <-> il. - pre_phase_factor = get_pre_phase_factor(k, q_vecs, atom_triplets); - real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[0], fc3, + real_to_reciprocal_elements(fc3_rec_elem, phase_factor1 + k * num_satom, + phase_factor0 + k * num_satom, fc3, is_compact_fc3, atom_triplets, k, j, i, 3); for (l = 0; l < 3; l++) { for (m = 0; m < 3; m++) { for (n = 0; n < 3; n++) { fc3_rec = phonoc_complex_prod( - fc3_rec_elem[n * 9 + m * 3 + l], pre_phase_factor); + fc3_rec_elem[n * 9 + m * 3 + l], pre_phase_factors[k]); fc3_reciprocal[(i * 3 + l) * num_band * num_band + (j * 3 + m) * num_band + k * 3 + n] = sum_lapack_complex_double( @@ -224,23 +273,19 @@ static void real_to_reciprocal_r0_average(lapack_complex_double *fc3_reciprocal, } } -static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem, - const double q1[3], const double q2[3], - const double *fc3, - const long is_compact_fc3, - const AtomTriplets *atom_triplets, - const long pi0, const long pi1, - const long pi2, const long leg_index) { +static void real_to_reciprocal_elements( + lapack_complex_double *fc3_rec_elem, + const lapack_complex_double *phase_factor1, + const lapack_complex_double *phase_factor2, const double *fc3, + const long is_compact_fc3, const AtomTriplets *atom_triplets, + const long pi0, const long pi1, 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; + long num_satom, adrs_shift; + lapack_complex_double phase_factor; double fc3_rec_real[27], fc3_rec_imag[27]; num_satom = atom_triplets->multi_dims[0]; - phase_factor2 = (lapack_complex_double *)malloc( - sizeof(lapack_complex_double) * num_satom); - for (i = 0; i < 27; i++) { fc3_rec_real[i] = 0; fc3_rec_imag[i] = 0; @@ -252,20 +297,11 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem, i = atom_triplets->p2s_map[pi0]; } - for (j = 0; j < num_satom; j++) { - adrs_vec2 = j * atom_triplets->multi_dims[1] + pi0; - phase_factor2[j] = get_phase_factor( - q2, atom_triplets->svecs, atom_triplets->multiplicity[adrs_vec2]); - } - for (j = 0; j < num_satom; j++) { if (atom_triplets->s2p_map[j] != atom_triplets->p2s_map[pi1]) { continue; } - adrs_vec1 = j * atom_triplets->multi_dims[1] + pi0; - phase_factor1 = get_phase_factor( - q1, atom_triplets->svecs, atom_triplets->multiplicity[adrs_vec1]); for (k = 0; k < num_satom; k++) { if (atom_triplets->s2p_map[k] != atom_triplets->p2s_map[pi2]) { continue; @@ -278,7 +314,8 @@ 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[k]); + phase_factor = + phonoc_complex_prod(phase_factor1[j], phase_factor2[k]); if ((leg_index == 1) && (atom_triplets->all_shortest[pi0 * num_satom * num_satom + @@ -308,9 +345,6 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem, fc3_rec_elem[i] = lapack_make_complex_double(fc3_rec_real[i], fc3_rec_imag[i]); } - - free(phase_factor2); - phase_factor2 = NULL; } // This function doesn't need to think about position + diff --git a/c/real_to_reciprocal.h b/c/real_to_reciprocal.h index d8bef345..3daa2807 100644 --- a/c/real_to_reciprocal.h +++ b/c/real_to_reciprocal.h @@ -52,5 +52,5 @@ void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal, const double q_vecs[3][3], const double *fc3, const long is_compact_fc3, const AtomTriplets *atom_triplets, - const long openmp_at_bands); + const long openmp_per_triplets); #endif diff --git a/c/reciprocal_to_normal.c b/c/reciprocal_to_normal.c index c3f39bef..bb2b360f 100644 --- a/c/reciprocal_to_normal.c +++ b/c/reciprocal_to_normal.c @@ -67,7 +67,7 @@ void reciprocal_to_normal_squared( const lapack_complex_double *eigvecs1, const lapack_complex_double *eigvecs2, const double *masses, const long *band_indices, const long num_band, - const double cutoff_frequency, const long openmp_at_bands) { + const double cutoff_frequency, const long openmp_per_triplets) { long i, j, num_atom; double *inv_sqrt_masses; lapack_complex_double *e0, *e1, *e2; @@ -95,7 +95,7 @@ void reciprocal_to_normal_squared( e2 = e1 + num_band * num_band; #ifdef _OPENMP -#pragma omp parallel for private(j) if (openmp_at_bands) +#pragma omp parallel for private(j) if (!openmp_per_triplets) #endif for (i = 0; i < num_band; i++) { for (j = 0; j < num_band; j++) { @@ -132,7 +132,7 @@ void reciprocal_to_normal_squared( #endif #ifdef _OPENMP -#pragma omp parallel for if (openmp_at_bands) +#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 && @@ -204,9 +204,8 @@ static double get_fc3_sum_blas(const lapack_complex_double *e0, const lapack_complex_double *e2, const lapack_complex_double *fc3_reciprocal, const long num_band) { - long i, j; + long i; lapack_complex_double *fc3_e12, *e_12, zero, one, retval; - const lapack_complex_double *fc3_i; e_12 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) * num_band * num_band); diff --git a/c/reciprocal_to_normal.h b/c/reciprocal_to_normal.h index 22bfa070..acdbdb30 100644 --- a/c/reciprocal_to_normal.h +++ b/c/reciprocal_to_normal.h @@ -45,6 +45,6 @@ void reciprocal_to_normal_squared( const lapack_complex_double *eigvecs1, const lapack_complex_double *eigvecs2, const double *masses, const long *band_indices, const long num_band, - const double cutoff_frequency, const long openmp_at_bands); + const double cutoff_frequency, const long openmp_per_triplets); #endif diff --git a/c/triplet.c b/c/triplet.c index a2adc238..27453fc2 100644 --- a/c/triplet.c +++ b/c/triplet.c @@ -64,7 +64,7 @@ void tpl_get_integration_weight( const long (*triplets)[3], const long num_triplets, const ConstBZGrid *bzgrid, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets, const long openmp_per_bands) { + const long openmp_per_triplets) { long i, num_band_prod; long tp_relative_grid_address[2][24][4][3]; @@ -82,7 +82,7 @@ void tpl_get_integration_weight( num_band0, tp_relative_grid_address, triplets[i], num_triplets, bzgrid, frequencies1, /* f1 */ num_band1, frequencies2, /* f2 */ - num_band2, tp_type, openmp_per_bands); + num_band2, tp_type, openmp_per_triplets); } } diff --git a/c/triplet.h b/c/triplet.h index b6489d9c..73c2ad51 100644 --- a/c/triplet.h +++ b/c/triplet.h @@ -66,7 +66,7 @@ void tpl_get_integration_weight( const long (*triplets)[3], const long num_triplets, const ConstBZGrid *bzgrid, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets, const long openmp_per_bands); + const long openmp_per_triplets); void tpl_get_integration_weight_with_sigma( double *iw, char *iw_zero, const double sigma, const double sigma_cutoff, const double *frequency_points, const long num_band0, diff --git a/c/triplet_iw.c b/c/triplet_iw.c index 4a88888a..b49b4e5d 100644 --- a/c/triplet_iw.c +++ b/c/triplet_iw.c @@ -67,7 +67,7 @@ void tpi_get_integration_weight( const long triplets[3], const long num_triplets, const ConstBZGrid *bzgrid, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_bands) { + const long openmp_per_triplets) { long max_i, j, b1, b2, b12, num_band_prod, adrs_shift; long vertices[2][24][4]; double g[3]; @@ -98,7 +98,7 @@ void tpi_get_integration_weight( #ifdef _OPENMP #pragma omp parallel for private(j, b1, b2, adrs_shift, g, \ - freq_vertices) if (openmp_per_bands) + freq_vertices) if (!openmp_per_triplets) #endif for (b12 = 0; b12 < num_band1 * num_band2; b12++) { b1 = b12 / num_band2; @@ -132,13 +132,13 @@ void tpi_get_integration_weight_with_sigma( double *iw, char *iw_zero, const double sigma, const double cutoff, const double *frequency_points, const long num_band0, const long triplet[3], const long const_adrs_shift, const double *frequencies, const long num_band, - const long tp_type, const long openmp_per_bands) { + const long tp_type, const long openmp_per_triplets) { long j, b12, b1, b2, adrs_shift; double f0, f1, f2, g0, g1, g2; #ifdef _OPENMP #pragma omp parallel for private(j, b1, b2, f0, f1, f2, g0, g1, g2, \ - adrs_shift) if (openmp_per_bands) + adrs_shift) if (!openmp_per_triplets) #endif for (b12 = 0; b12 < num_band * num_band; b12++) { b1 = b12 / num_band; diff --git a/c/triplet_iw.h b/c/triplet_iw.h index a2a6eacb..aacdae48 100644 --- a/c/triplet_iw.h +++ b/c/triplet_iw.h @@ -43,12 +43,12 @@ void tpi_get_integration_weight( const long triplets[3], const long num_triplets, const ConstBZGrid *bzgrid, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_bands); + const long openmp_per_triplets); void tpi_get_integration_weight_with_sigma( double *iw, char *iw_zero, const double sigma, const double cutoff, const double *frequency_points, const long num_band0, const long triplet[3], const long const_adrs_shift, const double *frequencies, const long num_band, - const long tp_type, const long openmp_per_bands); + const long tp_type, const long openmp_per_triplets); void tpi_get_neighboring_grid_points(long *neighboring_grid_points, const long grid_point, const long (*relative_grid_address)[3], From 8790ac4f02dd736bfacab6149825fc19e61184f8 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Thu, 18 Apr 2024 17:17:35 +0900 Subject: [PATCH 11/16] Optimize get_fc3_sum -> get_fc3_sum_blas_like --- c/imag_self_energy_with_g.c | 4 +- c/pp_collision.c | 5 +- c/reciprocal_to_normal.c | 121 ++++++++++++++++++++++++++++-------- 3 files changed, 98 insertions(+), 32 deletions(-) diff --git a/c/imag_self_energy_with_g.c b/c/imag_self_energy_with_g.c index d6d32e05..d9aac798 100644 --- a/c/imag_self_energy_with_g.c +++ b/c/imag_self_energy_with_g.c @@ -221,7 +221,7 @@ void ise_imag_self_energy_at_triplet( const long triplet[3], const long triplet_weight, const double *g1, const double *g2_3, const long (*g_pos)[4], const long num_g_pos, const double *temperatures, const long num_temps, - const double cutoff_frequency, const long openmp_possible, + const double cutoff_frequency, const long openmp_per_triplets, const long at_a_frequency_point) { long i, j; double *n1, *n2, *ise_at_g_pos; @@ -238,7 +238,7 @@ void ise_imag_self_energy_at_triplet( } #ifdef _OPENMP -#pragma omp parallel for private(j, g_pos_3) if (openmp_possible) +#pragma omp parallel for private(j, g_pos_3) if (!openmp_per_triplets) #endif for (i = 0; i < num_g_pos; i++) { if (at_a_frequency_point) { diff --git a/c/pp_collision.c b/c/pp_collision.c index 6258f083..9cd6b9e8 100644 --- a/c/pp_collision.c +++ b/c/pp_collision.c @@ -96,7 +96,6 @@ 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) @@ -234,12 +233,12 @@ static void get_collision( fc3_normal_squared, num_band0, num_band, g_pos, num_g_pos, frequencies, eigenvectors, triplet, bzgrid, fc3, is_compact_fc3, atom_triplets, masses, band_indices, symmetrize_fc3_q, cutoff_frequency, 0, 0, - 1 - openmp_per_triplets); + openmp_per_triplets); 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, 1 - openmp_per_triplets, 0); + num_temps, cutoff_frequency, openmp_per_triplets, 0); free(fc3_normal_squared); fc3_normal_squared = NULL; diff --git a/c/reciprocal_to_normal.c b/c/reciprocal_to_normal.c index bb2b360f..ca7a277c 100644 --- a/c/reciprocal_to_normal.c +++ b/c/reciprocal_to_normal.c @@ -59,6 +59,11 @@ static double get_fc3_sum_blas(const lapack_complex_double *e0, const lapack_complex_double *e2, const lapack_complex_double *fc3_reciprocal, const long num_band); +static double get_fc3_sum_blas_like(const lapack_complex_double *e0, + const lapack_complex_double *e1, + const lapack_complex_double *e2, + const lapack_complex_double *fc3_reciprocal, + const long num_band); void reciprocal_to_normal_squared( double *fc3_normal_squared, const long (*g_pos)[4], const long num_g_pos, const lapack_complex_double *fc3_reciprocal, const double *freqs0, @@ -68,7 +73,7 @@ void reciprocal_to_normal_squared( const lapack_complex_double *eigvecs2, const double *masses, const long *band_indices, const long num_band, const double cutoff_frequency, const long openmp_per_triplets) { - long i, j, num_atom; + long i, j, ij, num_atom, use_multithreaded_blas; double *inv_sqrt_masses; lapack_complex_double *e0, *e1, *e2; @@ -95,26 +100,26 @@ void reciprocal_to_normal_squared( e2 = e1 + num_band * num_band; #ifdef _OPENMP -#pragma omp parallel for private(j) if (!openmp_per_triplets) +#pragma omp parallel for private(i, j) if (!openmp_per_triplets) #endif - for (i = 0; i < num_band; i++) { - for (j = 0; j < num_band; j++) { - e0[i * num_band + j] = lapack_make_complex_double( - lapack_complex_double_real(eigvecs0[j * num_band + i]) * - inv_sqrt_masses[j], - lapack_complex_double_imag(eigvecs0[j * num_band + i]) * - inv_sqrt_masses[j]); - e1[i * num_band + j] = lapack_make_complex_double( - lapack_complex_double_real(eigvecs1[j * num_band + i]) * - inv_sqrt_masses[j], - lapack_complex_double_imag(eigvecs1[j * num_band + i]) * - inv_sqrt_masses[j]); - e2[i * num_band + j] = lapack_make_complex_double( - lapack_complex_double_real(eigvecs2[j * num_band + i]) * - inv_sqrt_masses[j], - lapack_complex_double_imag(eigvecs2[j * num_band + i]) * - inv_sqrt_masses[j]); - } + for (ij = 0; ij < num_band * num_band; ij++) { + i = ij / num_band; + j = ij % num_band; + e0[i * num_band + j] = lapack_make_complex_double( + lapack_complex_double_real(eigvecs0[j * num_band + i]) * + inv_sqrt_masses[j], + lapack_complex_double_imag(eigvecs0[j * num_band + i]) * + inv_sqrt_masses[j]); + e1[i * num_band + j] = lapack_make_complex_double( + lapack_complex_double_real(eigvecs1[j * num_band + i]) * + inv_sqrt_masses[j], + lapack_complex_double_imag(eigvecs1[j * num_band + i]) * + inv_sqrt_masses[j]); + e2[i * num_band + j] = lapack_make_complex_double( + lapack_complex_double_real(eigvecs2[j * num_band + i]) * + inv_sqrt_masses[j], + lapack_complex_double_imag(eigvecs2[j * num_band + i]) * + inv_sqrt_masses[j]); } free(inv_sqrt_masses); @@ -131,20 +136,39 @@ void reciprocal_to_normal_squared( loopStartCPUTime = clock(); #endif +#ifdef MULTITHREADED_BLAS + if (openmp_per_triplets) { + use_multithreaded_blas = 0; + } else { + use_multithreaded_blas = 1; + } +#else + use_multithreaded_blas = 0; #ifdef _OPENMP #pragma omp parallel for if (!openmp_per_triplets) +#endif #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(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]]); + if (use_multithreaded_blas) { + fc3_normal_squared[g_pos[i][3]] = + get_fc3_sum_blas(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]] = + 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; } @@ -234,3 +258,46 @@ static double get_fc3_sum_blas(const lapack_complex_double *e0, lapack_complex_double_imag(retval) * lapack_complex_double_imag(retval); } + +static double get_fc3_sum_blas_like(const lapack_complex_double *e0, + const lapack_complex_double *e1, + const lapack_complex_double *e2, + const lapack_complex_double *fc3_reciprocal, + const long num_band) { + long i, j; + double sum_real, sum_imag, retval_real, retval_imag; + lapack_complex_double *e_12, fc3_e_12, fc3_e_012; + + e_12 = (lapack_complex_double *)malloc(sizeof(lapack_complex_double) * + num_band * num_band); + + for (i = 0; i < num_band; i++) { + memcpy(e_12 + i * num_band, e2, 16 * num_band); + for (j = 0; j < num_band; j++) { + e_12[i * num_band + j] = + phonoc_complex_prod(e1[i], e_12[i * num_band + j]); + } + } + + retval_real = 0; + retval_imag = 0; + for (i = 0; i < num_band; i++) { + sum_real = 0; + sum_imag = 0; + for (j = 0; j < num_band * num_band; j++) { + fc3_e_12 = phonoc_complex_prod( + fc3_reciprocal[i * num_band * num_band + j], e_12[j]); + sum_real += lapack_complex_double_real(fc3_e_12); + sum_imag += lapack_complex_double_imag(fc3_e_12); + } + fc3_e_012 = phonoc_complex_prod( + e0[i], lapack_make_complex_double(sum_real, sum_imag)); + retval_real += lapack_complex_double_real(fc3_e_012); + retval_imag += lapack_complex_double_imag(fc3_e_012); + } + + free(e_12); + e_12 = NULL; + + return retval_real * retval_real + retval_imag * retval_imag; +} From 400fe29520a3d97384b11e6423a2322182b2829f Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Thu, 18 Apr 2024 18:04:18 +0900 Subject: [PATCH 12/16] Minor refactoring --- c/real_to_reciprocal.c | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/c/real_to_reciprocal.c b/c/real_to_reciprocal.c index 00f88068..29f9212b 100644 --- a/c/real_to_reciprocal.c +++ b/c/real_to_reciprocal.c @@ -149,15 +149,15 @@ static void real_to_reciprocal_legacy( const long openmp_per_triplets) { long i, j, k, l, m, n, ijk; long num_patom, num_satom, num_band; - lapack_complex_double fc3_rec_elem[27], fc3_rec; + lapack_complex_double fc3_rec_elem[27]; num_patom = atom_triplets->multi_dims[1]; num_satom = atom_triplets->multi_dims[0]; num_band = num_patom * 3; #ifdef _OPENMP -#pragma omp parallel for private(i, j, k, l, m, n, fc3_rec_elem, \ - fc3_rec) if (!openmp_per_triplets) +#pragma omp parallel for private(i, j, k, l, m, n, \ + fc3_rec_elem) if (!openmp_per_triplets) #endif for (ijk = 0; ijk < num_patom * num_patom * num_patom; ijk++) { i = ijk / (num_patom * num_patom); @@ -170,14 +170,10 @@ static void real_to_reciprocal_legacy( for (l = 0; l < 3; l++) { for (m = 0; m < 3; m++) { for (n = 0; n < 3; n++) { - fc3_rec = phonoc_complex_prod( - fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factors[i]); fc3_reciprocal[(i * 3 + l) * num_band * num_band + (j * 3 + m) * num_band + k * 3 + n] = - sum_lapack_complex_double( - fc3_reciprocal[(i * 3 + l) * num_band * num_band + - (j * 3 + m) * num_band + k * 3 + n], - fc3_rec); + phonoc_complex_prod(fc3_rec_elem[l * 9 + m * 3 + n], + pre_phase_factors[i]); } } } From ac874d115f0002fac149396266033785b22d2e2e Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Fri, 19 Apr 2024 14:47:15 +0900 Subject: [PATCH 13/16] Preparationg for v3 release --- doc/auxiliary-tools.md | 13 +- doc/changelog.md | 20 +- doc/conf.py | 4 +- phono3py/api_phono3py.py | 306 +------------------------ phono3py/cui/create_force_constants.py | 104 ++------- phono3py/cui/create_supercells.py | 25 -- phono3py/cui/kaccum_script.py | 7 +- phono3py/cui/load.py | 70 ++---- phono3py/cui/phono3py_script.py | 55 ++--- phono3py/file_IO.py | 133 +++-------- phono3py/other/isotope.py | 34 +-- phono3py/other/kaccum.py | 36 +-- phono3py/phonon3/dataset.py | 3 +- phono3py/phonon3/fc3.py | 2 +- phono3py/version.py | 2 +- test/other/test_kaccum.py | 93 ++------ 16 files changed, 150 insertions(+), 757 deletions(-) diff --git a/doc/auxiliary-tools.md b/doc/auxiliary-tools.md index 7432e6ec..9f6b3d7d 100644 --- a/doc/auxiliary-tools.md +++ b/doc/auxiliary-tools.md @@ -39,21 +39,18 @@ Let's computer lattice thermal conductivity of Si using the `Si-PBEsol` example found in the example directory. ```bash -% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="11 11 11" --sym-fc --br +% phono3py --mesh="11 11 11" --sym-fc --br ``` Then using the output file, `kappa-m111111.hdf5`, run `phono3py-kaccum` as follows: ```bash -% phono3py-kaccum --pa="F" -c POSCAR-unitcell kappa-m111111.hdf5 |tee kaccum.dat +% phono3py-kaccum kappa-m111111.hdf5 |tee kaccum.dat ``` -Here `--pa` is optional. The definition of `--pa` option is same as -{ref}`pa_option`. `POSCAR-unitcell` is the unit cell filename that is specified -with `-c` option. `kappa-m111111.hdf5` is the output file of thermal -conductivity calculation, which is passed to `phono3py-kaccum` as the first -argument. +`kappa-m111111.hdf5` is the output file of thermal conductivity calculation, +which is passed to `phono3py-kaccum` as the first argument. The format of the output is as follows: The first column gives frequency in THz, and the second to seventh columns give the cumulative lattice thermal @@ -99,7 +96,7 @@ Let `phono3py-kaccum` read a QE (pw) unit cell file with `-c` option, for example: ```bash -% phono3py-kaccum --qe --pa="F" -c Si.in kappa-m191919.hdf5 +% phono3py-kaccum --qe kappa-m191919.hdf5 ``` ```{image} Si-kaccum-pwscf.png diff --git a/doc/changelog.md b/doc/changelog.md index c6911716..e3fb9719 100644 --- a/doc/changelog.md +++ b/doc/changelog.md @@ -2,6 +2,24 @@ # Change Log +## Apr-19-2024: Version 3.0.0 + +This is a major version release. There are backward-incompatible changes. + +- Method to transform supercell third-order force constants fc3 in real to + reciprocal space was changed as described at Version 2.9.0 changelog below. + This results in the change of results with respect to those obtained by + phono3py version 2. To emulate v2 behaviour, use `--v2` option in phono3py + command line script. For `Phono3py` class , `make_r0_average=True` (default) + when instantiating it, and similarly for `phono3py.load` function. +- Completely dropped support of `disp_fc3.yaml` and `disp_fc2.yaml`. +- Dropped support of old style usage of `phono3py-kaccum` script. +- Removed functions `write_fc3_dat`, `write_triplets`, `write_grid_address` in + `file_IO.py`. +- Removed methods marked by `DeprecationWarning`. +- Removed `masses`, `band_indices`, `sigmas`, `sigma_cutoff` parameters from + `Phono3py.__init__()`. + ## Mar-20-2024: Version 2.10.0 - Maintenance release @@ -100,7 +118,7 @@ ## Jul-22-2021: Version 2.0.0 -This is a major version release. There are some backword-incompatible changes. +This is a major version release. There are some backward-incompatible changes. 1. Grid point indexing system to address grid points of q-points is changed. 2. Array data types of most of the integer arrays are changed to `dtype='int_'` diff --git a/doc/conf.py b/doc/conf.py index 140fdb53..c0730891 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -58,9 +58,9 @@ # built documents. # # The short X.Y version. -version = "2.10" +version = "3.0" # The full version, including alpha/beta/rc tags. -release = "2.10.0" +release = "3.0.0" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index 86efd8a7..8de413d4 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -34,7 +34,6 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -import warnings from collections.abc import Sequence from typing import Optional, Union @@ -144,17 +143,13 @@ def __init__( supercell_matrix=None, primitive_matrix=None, phonon_supercell_matrix=None, - masses=None, - band_indices=None, - sigmas=None, - sigma_cutoff=None, cutoff_frequency=1e-4, frequency_factor_to_THz=VaspToTHz, is_symmetry=True, is_mesh_symmetry=True, use_grg=False, SNF_coordinates="reciprocal", - make_r0_average: bool = False, + make_r0_average: bool = True, symprec=1e-5, calculator: Optional[str] = None, log_level=0, @@ -185,14 +180,6 @@ def __init__( than that of fc3. Unless setting this, supercell_matrix is used. This is only valide when unitcell or unitcell_filename is given. Default is None. - masses : Deprecated. - Use Phono3py.masses attribute after instanciation. - band_indices : Deprecated. - Use Phono3py.band_indices attribute after instanciation. - sigmas : Deprecated. - Use Phono3py.sigmas attribute after instanciation. - sigma_cutoff : Deprecated. - Use Phono3py.sigma_cutoff attribute after instanciation. cutoff_frequency : float, optional Phonon frequency below this value is ignored when the cutoff is needed for the computation. Default is 1e-4. @@ -308,38 +295,6 @@ def __init__( self._band_indices_flatten = None self._set_band_indices() - if masses is not None: - warnings.warn( - "Phono3py init parameter of masses is deprecated. " - "Use Phono3py.masses attribute instead.", - DeprecationWarning, - ) - self.masses = masses - - if band_indices is not None: - warnings.warn( - "Phono3py init parameter of band_indices is deprecated. " - "Use Phono3py.band_indices attribute instead.", - DeprecationWarning, - ) - self.band_indices = band_indices - - if sigmas is not None: - warnings.warn( - "Phono3py init parameter of sigmas is deprecated. " - "Use Phono3py.sigmas attribute instead.", - DeprecationWarning, - ) - self.sigmas = sigmas - - if sigma_cutoff is not None: - warnings.warn( - "Phono3py init parameter of sigma_cutoff is deprecated. " - "Use Phono3py.sigma_cutoff attribute instead.", - DeprecationWarning, - ) - self.sigma_cutoff = sigma_cutoff - @property def version(self): """Return phono3py release version number. @@ -350,15 +305,6 @@ def version(self): """ return __version__ - def get_version(self): - """Return phono3py release version number.""" - warnings.warn( - "Phono3py.get_version() is deprecated." - "Use Phono3py.version attribute instead.", - DeprecationWarning, - ) - return self.version - @property def calculator(self) -> Optional[str]: """Return calculator interface name. @@ -386,22 +332,6 @@ def fc3(self): def fc3(self, fc3): self._fc3 = fc3 - def get_fc3(self): - """Return third order force constants (fc3).""" - warnings.warn( - "Phono3py.get_fc3() is deprecated." "Use Phono3py.fc3 attribute instead.", - DeprecationWarning, - ) - return self.fc3 - - def set_fc3(self, fc3): - """Set fc3.""" - warnings.warn( - "Phono3py.set_fc3() is deprecated." "Use Phono3py.fc3 attribute instead.", - DeprecationWarning, - ) - self.fc3 = fc3 - @property def fc2(self): """Setter and getter of second order force constants (fc2). @@ -419,22 +349,6 @@ def fc2(self): def fc2(self, fc2): self._fc2 = fc2 - def get_fc2(self): - """Return second order force constants (fc2).""" - warnings.warn( - "Phono3py.get_fc2() is deprecated." "Use Phono3py.fc2 attribute instead.", - DeprecationWarning, - ) - return self.fc2 - - def set_fc2(self, fc2): - """Set fc2.""" - warnings.warn( - "Phono3py.set_fc2() is deprecated." "Use Phono3py.fc2 attribute instead.", - DeprecationWarning, - ) - self.fc2 = fc2 - @property def force_constants(self): """Return fc2. This is same as the getter attribute `fc2`.""" @@ -511,24 +425,6 @@ def nac_params(self, nac_params): if self._interaction is not None: self._init_dynamical_matrix() - def get_nac_params(self): - """Return NAC parameters.""" - warnings.warn( - "Phono3py.get_nac_params() is deprecated." - "Use Phono3py.nac_params attribute instead.", - DeprecationWarning, - ) - return self.nac_params - - def set_nac_params(self, nac_params): - """Set NAC parameters.""" - warnings.warn( - "Phono3py.set_nac_params() is deprecated." - "Use Phono3py.nac_params attribute instead.", - DeprecationWarning, - ) - self.nac_params = nac_params - @property def dynamical_matrix(self): """Return DynamicalMatrix instance. @@ -552,15 +448,6 @@ def primitive(self) -> Primitive: """ return self._primitive - def get_primitive(self): - """Return primitive cell.""" - warnings.warn( - "Phono3py.get_primitive() is deprecated." - "Use Phono3py.primitive attribute instead.", - DeprecationWarning, - ) - return self.primitive - @property def unitcell(self) -> PhonopyAtoms: """Return Unit cell. @@ -571,15 +458,6 @@ def unitcell(self) -> PhonopyAtoms: """ return self._unitcell - def get_unitcell(self): - """Return Unit cell.""" - warnings.warn( - "Phono3py.get_unitcell() is deprecated." - "Use Phono3py.unitcell attribute instead.", - DeprecationWarning, - ) - return self.unitcell - @property def supercell(self) -> Supercell: """Return supercell. @@ -590,15 +468,6 @@ def supercell(self) -> Supercell: """ return self._supercell - def get_supercell(self): - """Return supercell.""" - warnings.warn( - "Phono3py.get_supercell() is deprecated." - "Use Phono3py.supercell attribute instead.", - DeprecationWarning, - ) - return self.supercell - @property def phonon_supercell(self) -> Supercell: """Return supercell for fc2. @@ -609,15 +478,6 @@ def phonon_supercell(self) -> Supercell: """ return self._phonon_supercell - def get_phonon_supercell(self): - """Return supercell for fc2.""" - warnings.warn( - "Phono3py.get_phonon_supercell() is deprecated." - "Use Phono3py.phonon_supercell attribute instead.", - DeprecationWarning, - ) - return self.phonon_supercell - @property def phonon_primitive(self) -> Primitive: """Return primitive cell for fc2. @@ -630,15 +490,6 @@ def phonon_primitive(self) -> Primitive: """ return self._phonon_primitive - def get_phonon_primitive(self): - """Return primitive cell for fc2.""" - warnings.warn( - "Phono3py.get_phonon_primitive() is deprecated." - "Use Phono3py.phonon_primitive attribute instead.", - DeprecationWarning, - ) - return self.phonon_primitive - @property def symmetry(self) -> Symmetry: """Return symmetry of supercell. @@ -649,15 +500,6 @@ def symmetry(self) -> Symmetry: """ return self._symmetry - def get_symmetry(self): - """Return symmetry of supercell.""" - warnings.warn( - "Phono3py.get_symmetry() is deprecated." - "Use Phono3py.symmetry attribute instead.", - DeprecationWarning, - ) - return self.symmetry - @property def primitive_symmetry(self) -> Symmetry: """Return symmetry of primitive cell. @@ -668,15 +510,6 @@ def primitive_symmetry(self) -> Symmetry: """ return self._primitive_symmetry - def get_primitive_symmetry(self): - """Return symmetry of primitive cell.""" - warnings.warn( - "Phono3py.get_primitive_symmetry() is deprecated." - "Use Phono3py.primitive_symmetry attribute instead.", - DeprecationWarning, - ) - return self.primitive_symmetry - @property def phonon_supercell_symmetry(self) -> Symmetry: """Return symmetry of supercell for fc2. @@ -687,15 +520,6 @@ def phonon_supercell_symmetry(self) -> Symmetry: """ return self._phonon_supercell_symmetry - def get_phonon_supercell_symmetry(self): - """Return symmetry of supercell for fc2.""" - warnings.warn( - "Phono3py.get_phonon_supercell_symmetry() is deprecated." - "Use Phono3py.phonon_supercell_symmetry attribute instead.", - DeprecationWarning, - ) - return self.phonon_supercell_symmetry - @property def supercell_matrix(self): """Return transformation matrix to supercell cell from unit cell. @@ -707,15 +531,6 @@ def supercell_matrix(self): """ return self._supercell_matrix - def get_supercell_matrix(self): - """Return transformation matrix to supercell cell from unit cell.""" - warnings.warn( - "Phono3py.get_supercell_matrix() is deprecated." - "Use Phono3py.supercell_matrix attribute instead.", - DeprecationWarning, - ) - return self.supercell_matrix - @property def phonon_supercell_matrix(self): """Return transformation matrix to phonon supercell from unit cell. @@ -727,15 +542,6 @@ def phonon_supercell_matrix(self): """ return self._phonon_supercell_matrix - def get_phonon_supercell_matrix(self): - """Return transformation matrix to phonon supercell from unit cell.""" - warnings.warn( - "Phono3py.get_phonon_supercell_matrix() is deprecated." - "Use Phono3py.phonon_supercell_matrix attribute instead.", - DeprecationWarning, - ) - return self.phonon_supercell_matrix - @property def primitive_matrix(self): """Return transformation matrix to primitive cell from unit cell. @@ -747,15 +553,6 @@ def primitive_matrix(self): """ return self._primitive_matrix - def get_primitive_matrix(self): - """Return transformation matrix to primitive cell from unit cell.""" - warnings.warn( - "Phono3py.get_primitive_matrix() is deprecated." - "Use Phono3py.primitive_matrix attribute instead.", - DeprecationWarning, - ) - return self.primitive_matrix - @property def unit_conversion_factor(self): """Return phonon frequency unit conversion factor. @@ -768,15 +565,6 @@ def unit_conversion_factor(self): """ return self._frequency_factor_to_THz - def set_displacement_dataset(self, dataset): - """Set displacement-force dataset.""" - warnings.warn( - "Phono3py.set_displacement_dataset() is deprecated." - "Use Phono3py.dataset attribute instead.", - DeprecationWarning, - ) - self._dataset = dataset - @property def dataset(self): """Setter and getter of displacement-force dataset. @@ -823,24 +611,6 @@ def dataset(self, dataset): self._supercells_with_displacements = None self._phonon_supercells_with_displacements = None - @property - def displacement_dataset(self): - """Return displacement-force dataset.""" - warnings.warn( - "Phono3py.displacement_dataset is deprecated." "Use Phono3py.dataset.", - DeprecationWarning, - ) - return self.dataset - - def get_displacement_dataset(self): - """Return displacement-force dataset.""" - warnings.warn( - "Phono3py.get_displacement_dataset() is deprecated." - "Use Phono3py.dataset.", - DeprecationWarning, - ) - return self.displacement_dataset - @property def phonon_dataset(self): """Setter and getter of displacement-force dataset for fc2. @@ -869,25 +639,6 @@ def phonon_dataset(self): def phonon_dataset(self, dataset): self._phonon_dataset = dataset - @property - def phonon_displacement_dataset(self): - """Return phonon dispalcement-force dataset.""" - warnings.warn( - "Phono3py.phonon_displacement_dataset is deprecated." - "Use Phono3py.phonon_dataset.", - DeprecationWarning, - ) - return self._phonon_dataset - - def get_phonon_displacement_dataset(self): - """Return phonon dispalcement-force dataset.""" - warnings.warn( - "Phono3py.get_phonon_displacement_dataset() is deprecated." - "Use Phono3py.phonon_dataset.", - DeprecationWarning, - ) - return self.phonon_displacement_dataset - @property def band_indices(self): """Setter and getter of band indices. @@ -903,15 +654,6 @@ def band_indices(self): def band_indices(self, band_indices): self._set_band_indices(band_indices=band_indices) - def set_band_indices(self, band_indices): - """Set band indices.""" - warnings.warn( - "Phono3py.set_band_indices() is deprecated." - "Use Phono3py.band_indices attribute instead.", - DeprecationWarning, - ) - self.band_indices = band_indices - def _set_band_indices(self, band_indices=None): if band_indices is None: num_band = len(self._primitive) * 3 @@ -955,15 +697,6 @@ def supercells_with_displacements(self): self._build_supercells_with_displacements() return self._supercells_with_displacements - def get_supercells_with_displacements(self): - """Return supercells with displacements.""" - warnings.warn( - "Phono3py.get_supercells_with_displacements() is deprecated." - "Use Phono3py.supercells_with_displacements attribute instead.", - DeprecationWarning, - ) - return self.supercells_with_displacements - @property def phonon_supercells_with_displacements(self): """Return supercells with displacements for fc2. @@ -982,16 +715,6 @@ def phonon_supercells_with_displacements(self): ) return self._phonon_supercells_with_displacements - def get_phonon_supercells_with_displacements(self): - """Return supercells with displacements for fc2.""" - warnings.warn( - "Phono3py.get_phonon_supercells_with_displacements() " - "is deprecated. Use Phono3py.phonon_supercells_with_displacements " - "attribute instead.", - DeprecationWarning, - ) - return self.phonon_supercells_with_displacements - @property def mesh_numbers(self): """Setter and getter of sampling mesh numbers in reciprocal space.""" @@ -1009,15 +732,6 @@ def thermal_conductivity(self): """Return thermal conductivity class instance.""" return self._thermal_conductivity - def get_thermal_conductivity(self): - """Return thermal conductivity class instance.""" - warnings.warn( - "Phono3py.get_thermal_conductivity() is deprecated." - "Use Phono3py.thermal_conductivity attribute instead.", - DeprecationWarning, - ) - return self.thermal_conductivity - @property def displacements(self): """Setter and getter displacements in supercells. @@ -1261,15 +975,6 @@ def phph_interaction(self): """Return Interaction instance.""" return self._interaction - def get_phph_interaction(self): - """Return Interaction instance.""" - warnings.warn( - "Phono3py.get_phph_interaction() is deprecated." - "Use Phono3py.phph_interaction attribute instead.", - DeprecationWarning, - ) - return self.phph_interaction - @property def detailed_gammas(self): """Return detailed gamma.""" @@ -1880,15 +1585,6 @@ def run_imag_self_energy( return vals - def write_imag_self_energy(self, filename=None): - """Write imaginary part of self energy to a file.""" - warnings.warn( - "Phono3py.write_imag_self_energy is deprecated." - "Use Phono3py.run_imag_self_energy with write_txt=True.", - DeprecationWarning, - ) - self._write_imag_self_energy(output_filename=filename) - def _write_imag_self_energy(self, output_filename=None): write_imag_self_energy( self._gammas, diff --git a/phono3py/cui/create_force_constants.py b/phono3py/cui/create_force_constants.py index d0042586..6165d405 100644 --- a/phono3py/cui/create_force_constants.py +++ b/phono3py/cui/create_force_constants.py @@ -54,8 +54,6 @@ from phono3py.cui.show_log import show_phono3py_force_constants_settings from phono3py.file_IO import ( get_length_of_first_line, - parse_disp_fc2_yaml, - parse_disp_fc3_yaml, parse_FORCES_FC2, parse_FORCES_FC3, read_fc2_from_hdf5, @@ -117,7 +115,6 @@ def create_phono3py_force_constants( ph3py_yaml, phono3py_yaml_filename, symmetrize_fc3r, - input_filename, settings.is_compact_fc, settings.cutoff_pair_distance, fc_calculator, @@ -167,7 +164,6 @@ def create_phono3py_force_constants( phono3py, ph3py_yaml, symmetrize_fc2, - input_filename, settings.is_compact_fc, fc_calculator, fc_calculator_options, @@ -178,7 +174,6 @@ def create_phono3py_force_constants( phono3py, ph3py_yaml, symmetrize_fc2, - input_filename, settings.is_compact_fc, fc_calculator, fc_calculator_options, @@ -207,7 +202,6 @@ def parse_forces( ph3py_yaml: Optional[Phono3pyYaml] = None, cutoff_pair_distance=None, force_filename: str = "FORCES_FC3", - disp_filename: Optional[str] = None, phono3py_yaml_filename: Optional[str] = None, fc_type=None, log_level=0, @@ -269,32 +263,7 @@ def parse_forces( force_to_eVperA=physical_units["force_to_eVperA"], ) - # No dataset is found. "disp_fc*.yaml" is read. But this is deprecated. - if dataset is None: - if disp_filename is None: - msg = ( - "Displacement dataset corresponding to " - f'"{force_filename}" not found.' - ) - raise RuntimeError(msg) - # Displacement dataset is obtained from disp_filename. - # can emit FileNotFoundError. - dataset = _read_disp_fc_yaml(disp_filename, fc_type) - filename_read_from = disp_filename - - # No forces should exist. Therefore only unit of displacements is - # converted. - if calculator is None: - if log_level: - print( - f'Displacements are read from "{disp_filename}", but ' - " the unit has not converted." - ) - else: - _convert_unit_in_dataset( - dataset, - distance_to_A=physical_units["distance_to_A"], - ) + assert dataset is not None if "natom" in dataset and dataset["natom"] != natom: msg = ( @@ -369,15 +338,6 @@ def get_fc_calculator_params(settings): return fc_calculator, fc_calculator_options -def _read_disp_fc_yaml(disp_filename, fc_type): - if fc_type == "phonon_fc2": - dataset = parse_disp_fc2_yaml(filename=disp_filename) - else: - dataset = parse_disp_fc3_yaml(filename=disp_filename) - - return dataset - - def _read_phono3py_fc3(phono3py, symmetrize_fc3r, input_filename, log_level): if input_filename is None: filename = "fc3.hdf5" @@ -469,13 +429,12 @@ def _create_phono3py_fc3( phono3py: Phono3py, ph3py_yaml: Optional[Phono3pyYaml], phono3py_yaml_filename: Optional[str], - symmetrize_fc3r, - input_filename, - is_compact_fc, - cutoff_pair_distance, - fc_calculator, - fc_calculator_options, - log_level, + symmetrize_fc3r: bool, + is_compact_fc: bool, + cutoff_pair_distance: Optional[float], + fc_calculator: Optional[str], + fc_calculator_options: Optional[str], + log_level: int, ): """Read or calculate fc3. @@ -491,14 +450,7 @@ def _create_phono3py_fc3( when the former value is smaller than the later. """ - # disp_fc3.yaml is obsolete. - if input_filename is None: - disp_filename = "disp_fc3.yaml" - else: - disp_filename = "disp_fc3." + input_filename + ".yaml" - - # If disp_fc3.yaml is not found, default phono3py.yaml file is - _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml) + _ph3py_yaml = _get_default_ph3py_yaml(ph3py_yaml) try: dataset = parse_forces( @@ -506,7 +458,6 @@ def _create_phono3py_fc3( ph3py_yaml=_ph3py_yaml, cutoff_pair_distance=cutoff_pair_distance, force_filename="FORCES_FC3", - disp_filename=disp_filename, phono3py_yaml_filename=phono3py_yaml_filename, fc_type="fc3", log_level=log_level, @@ -534,25 +485,18 @@ def _create_phono3py_fc2( phono3py: Phono3py, ph3py_yaml: Optional[Phono3pyYaml], symmetrize_fc2, - input_filename, is_compact_fc, fc_calculator, fc_calculator_options, log_level, ): - if input_filename is None: - disp_filename = "disp_fc3.yaml" - else: - disp_filename = "disp_fc3." + input_filename + ".yaml" - - _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml) + _ph3py_yaml = _get_default_ph3py_yaml(ph3py_yaml) try: dataset = parse_forces( phono3py, ph3py_yaml=_ph3py_yaml, force_filename="FORCES_FC3", - disp_filename=disp_filename, fc_type="fc2", log_level=log_level, ) @@ -573,40 +517,30 @@ def _create_phono3py_fc2( ) -def _get_ph3py_yaml(disp_filename, ph3py_yaml: Optional[Phono3pyYaml]): +def _get_default_ph3py_yaml(ph3py_yaml: Optional[Phono3pyYaml]): _ph3py_yaml = ph3py_yaml - # Try to use phono3py.phonon_dataset when the disp file not found - if not os.path.isfile(disp_filename): - disp_filename = None - if _ph3py_yaml is None and os.path.isfile("phono3py_disp.yaml"): - _ph3py_yaml = Phono3pyYaml() - _ph3py_yaml.read("phono3py_disp.yaml") + if _ph3py_yaml is None and os.path.isfile("phono3py_disp.yaml"): + _ph3py_yaml = Phono3pyYaml() + _ph3py_yaml.read("phono3py_disp.yaml") return _ph3py_yaml def _create_phono3py_phonon_fc2( phono3py: Phono3py, ph3py_yaml: Optional[Phono3pyYaml], - symmetrize_fc2, - input_filename, - is_compact_fc, - fc_calculator, - fc_calculator_options, - log_level, + symmetrize_fc2: bool, + is_compact_fc: bool, + fc_calculator: Optional[str], + fc_calculator_options: Optional[str], + log_level: int, ): - if input_filename is None: - disp_filename = "disp_fc2.yaml" - else: - disp_filename = "disp_fc2." + input_filename + ".yaml" - - _ph3py_yaml = _get_ph3py_yaml(disp_filename, ph3py_yaml) + _ph3py_yaml = _get_default_ph3py_yaml(ph3py_yaml) try: dataset = parse_forces( phono3py, ph3py_yaml=_ph3py_yaml, force_filename="FORCES_FC2", - disp_filename=disp_filename, fc_type="phonon_fc2", log_level=log_level, ) diff --git a/phono3py/cui/create_supercells.py b/phono3py/cui/create_supercells.py index 55977358..1b6cbbc4 100644 --- a/phono3py/cui/create_supercells.py +++ b/phono3py/cui/create_supercells.py @@ -37,7 +37,6 @@ from phonopy.interface.calculator import write_supercells_with_displacements from phono3py import Phono3py -from phono3py.file_IO import write_disp_fc2_yaml, write_disp_fc3_yaml from phono3py.interface.calculator import ( get_additional_info_to_write_fc2_supercells, get_additional_info_to_write_supercells, @@ -52,18 +51,12 @@ def create_phono3py_supercells( output_filename=None, interface_mode="vasp", log_level=1, - write_disp_yaml=False, ): """Create displacements and supercells. Distance unit used is that for the calculator interface. The default unit is Angstron. - Parameters - ---------- - write_disp_yaml : bool - Write old-style files of disp_fc3.yaml and disp_fc2.yaml. Default is False. - """ optional_structure_info = cell_info["optional_structure_info"] @@ -92,15 +85,6 @@ def create_phono3py_supercells( print('Unit cell was read from "%s".' % optional_structure_info[0]) print("Displacement distance: %s" % distance) - if write_disp_yaml: - if output_filename is None: - filename = "disp_fc3.yaml" - else: - filename = "disp_fc3." + output_filename + ".yaml" - num_disps, num_disp_files = write_disp_fc3_yaml( - phono3py.dataset, phono3py.supercell, filename=filename - ) - ids = [] disp_cells = [] for i, cell in enumerate(phono3py.supercells_with_displacements): @@ -132,15 +116,6 @@ def create_phono3py_supercells( print("Number of displacement supercell files created: %d" % num_disp_files) if phono3py.phonon_supercell_matrix is not None: - if write_disp_yaml: - if output_filename is None: - filename = "disp_fc2.yaml" - else: - filename = "disp_fc2." + output_filename + ".yaml" - num_disps = write_disp_fc2_yaml( - phono3py.phonon_dataset, phono3py.phonon_supercell, filename=filename - ) - additional_info = get_additional_info_to_write_fc2_supercells( interface_mode, phono3py.phonon_supercell_matrix ) diff --git a/phono3py/cui/kaccum_script.py b/phono3py/cui/kaccum_script.py index 744f8bbe..0937c3e7 100644 --- a/phono3py/cui/kaccum_script.py +++ b/phono3py/cui/kaccum_script.py @@ -358,9 +358,10 @@ def main(): """ args = _get_parser() primitive = None - if len(args.filenames) > 1: # deprecated - cell, f_kappa = _read_files(args) - primitive_matrix = _analyze_primitive_matrix_option(args, unitcell=cell) + if len(args.filenames) > 1: + raise RuntimeError( + 'Use of "phono3py-kaccum CRYSTAL_STRUCTURE_FILE" is not supported.' + ) else: interface_mode = _get_calculator(args) cell_info = _read_files_by_collect_cell_info(args.cell_filename, interface_mode) diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py index 59456ed6..d9ae3c8a 100644 --- a/phono3py/cui/load.py +++ b/phono3py/cui/load.py @@ -94,17 +94,17 @@ def load( can be overwritten. 'fc3.hdf5' is read if found in current directory. Unless 'fc3.hdf5' is found - and if 'FORCES_FC3' and 'disp_fc3.yaml" are found, these are read and fc3 - and fc2 are produced. + and if 'FORCES_FC3' and 'phono3py_disp.yaml" are found, these are read and + fc3 and fc2 are produced. if 'fc2.hdf5' is found, this is read. Unless 'fc2.hdf5' is found and if - 'FORCES_FC2' and 'disp_fc2.yaml" are found, these are read and fc2 is + 'FORCES_FC2' and 'phono3py_disp.yaml" are found, these are read and fc2 is produced. When force_sets_filename and force_constants_filename are not given, 'FORCES_FC3' and 'FORCES_FC2' are looked for in the current directory as the default behaviour. When 'FORCES_FC3' ('FORCES_FC2') is given in the type-1 - format, 'disp_fc3.yaml' ('disp_fc2.yaml') is also necessary and read. + format, 'phono3py_disp.yaml' is also necessary and read. Crystal structure ----------------- @@ -121,12 +121,12 @@ def load( priority: 1. fc3_filename (fc2_filename) 2. forces_fc3_filename (forces_fc2_filename). Do not forget that for - type-1 format, disp_fc3.yaml (disp_fc2.yaml) has to be given, too. + type-1 format, phono3py_disp.yaml has to be given, too. 3. 'fc3.hdf5' and 'fc2.hdf5' are searched in current directory. 4. 'FORCES_FC3' and 'FORCES_FC2' are searched in current directory. - 'FORCES_FC2' is optional. For type-1 format, 'disp_fc3.yaml' and - optionally 'disp_fc2.yaml' are also searched in current directory. - When 'FORCES_FC2' is not found, 'FORCES_FC3' is used to create fc2. + 'FORCES_FC2' is optional. For type-1 format, 'phono3py_disp.yaml' is + also searched in current directory. When 'FORCES_FC2' is not found, + 'FORCES_FC3' is used to create fc2. Parameters for non-analytical term correctiion (NAC) ---------------------------------------------------- @@ -188,11 +188,11 @@ def load( correction parameters. forces_fc3_filename : sequence or os.PathLike, optional A two-elemental sequence of filenames corresponding to ('FORCES_FC3', - 'disp_fc3.yaml') in the type-1 format or a filename (os.PathLike) + 'phono3py_disp.yaml') in the type-1 format or a filename (os.PathLike) corresponding to 'FORCES_FC3' in the type-2 format. Default is None. forces_fc2_filename : os.PathLike or sequence, optional A two-elemental sequence of filenames corresponding to ('FORCES_FC2', - 'disp_fc2.yaml') in the type-1 format or a filename (os.PathLike) + 'phono3py_disp.yaml') in the type-1 format or a filename (os.PathLike) corresponding to 'FORCES_FC2' in the type-2 format. Default is None. fc3_filename : os.PathLike, optional Filename of a file corresponding to 'fc3.hdf5', a file contains @@ -233,10 +233,10 @@ def load( use_grg : bool, optional Use generalized regular grid when True. Default is False. make_r0_average : bool, optional - fc3 transformation from real to reciprocal space is done - around three atoms and averaged when True. Default is False, i.e., - only around the first atom. Setting False is for rough compatibility - with v2.x. Default is True. + fc3 transformation from real to reciprocal space is done around three + atoms and averaged when True. Default is False, i.e., only around the + first atom. Setting False is for rough compatibility with v2.x. Default + is True. symprec : float, optional Tolerance used to find crystal symmetry. Default is 1e-5. log_level : int, optional @@ -477,16 +477,11 @@ def _set_dataset_or_fc3( if log_level: print('fc3 was read from "%s".' % fc3_filename) elif forces_fc3_filename is not None: - if isinstance(forces_fc3_filename, os.PathLike): - force_filename = forces_fc3_filename - disp_filename = None - else: - force_filename, disp_filename = forces_fc3_filename + force_filename = forces_fc3_filename _set_dataset_for_fc3( ph3py, ph3py_yaml, force_filename, - disp_filename, phono3py_yaml_filename, cutoff_pair_distance, log_level, @@ -499,17 +494,10 @@ def _set_dataset_or_fc3( if log_level: print('fc3 was read from "fc3.hdf5".') elif os.path.isfile("FORCES_FC3"): - disp_filename = None - if os.path.isfile("disp_fc3.yaml"): - if ph3py_yaml is None: - disp_filename = "disp_fc3.yaml" - elif ph3py_yaml.dataset is None: - disp_filename = "disp_fc3.yaml" _set_dataset_for_fc3( ph3py, ph3py_yaml, "FORCES_FC3", - disp_filename, phono3py_yaml_filename, cutoff_pair_distance, log_level, @@ -523,7 +511,6 @@ def _set_dataset_or_fc3( ph3py, ph3py_yaml, None, - None, phono3py_yaml_filename, cutoff_pair_distance, log_level, @@ -548,16 +535,11 @@ def _set_dataset_phonon_dataset_or_fc2( if log_level: print('fc2 was read from "%s".' % fc2_filename) elif forces_fc2_filename is not None: - if isinstance(forces_fc2_filename, os.PathLike): - force_filename = forces_fc2_filename - disp_filename = None - else: - force_filename, disp_filename = forces_fc2_filename + force_filename = forces_fc2_filename _set_dataset_for_fc2( ph3py, ph3py_yaml, force_filename, - disp_filename, "phonon_fc2", log_level, ) @@ -569,17 +551,10 @@ def _set_dataset_phonon_dataset_or_fc2( if log_level: print('fc2 was read from "fc2.hdf5".') elif os.path.isfile("FORCES_FC2"): - disp_filename = None - if os.path.isfile("disp_fc2.yaml"): - if ph3py_yaml is None: - disp_filename = "disp_fc2.yaml" - elif ph3py_yaml.phonon_dataset is None: - disp_filename = "disp_fc2.yaml" _set_dataset_for_fc2( ph3py, ph3py_yaml, "FORCES_FC2", - disp_filename, "phonon_fc2", log_level, ) @@ -592,7 +567,6 @@ def _set_dataset_phonon_dataset_or_fc2( ph3py, ph3py_yaml, None, - None, "phonon_fc2", log_level, ) @@ -605,23 +579,15 @@ def _set_dataset_phonon_dataset_or_fc2( ph3py, ph3py_yaml, None, - None, "fc2", log_level, ) elif os.path.isfile("FORCES_FC3"): # suppose fc3.hdf5 is read but fc2.hdf5 doesn't exist. - disp_filename = None - if os.path.isfile("disp_fc3.yaml"): - if ph3py_yaml is None: - disp_filename = "disp_fc3.yaml" - elif ph3py_yaml.dataset is None: - disp_filename = "disp_fc3.yaml" _set_dataset_for_fc2( ph3py, ph3py_yaml, "FORCES_FC3", - disp_filename, "fc2", log_level, ) @@ -632,7 +598,6 @@ def _set_dataset_for_fc3( ph3py: Phono3py, ph3py_yaml: Optional[Phono3pyYaml], force_filename, - disp_filename, phono3py_yaml_filename, cutoff_pair_distance, log_level, @@ -642,7 +607,6 @@ def _set_dataset_for_fc3( ph3py_yaml=ph3py_yaml, cutoff_pair_distance=cutoff_pair_distance, force_filename=force_filename, - disp_filename=disp_filename, phono3py_yaml_filename=phono3py_yaml_filename, fc_type="fc3", log_level=log_level, @@ -653,7 +617,6 @@ def _set_dataset_for_fc2( ph3py: Phono3py, ph3py_yaml: Optional[Phono3pyYaml], force_filename, - disp_filename, fc_type, log_level, ): @@ -661,7 +624,6 @@ def _set_dataset_for_fc2( ph3py, ph3py_yaml=ph3py_yaml, force_filename=force_filename, - disp_filename=disp_filename, fc_type=fc_type, log_level=log_level, ) diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py index 86ffeba2..138e6e5d 100644 --- a/phono3py/cui/phono3py_script.py +++ b/phono3py/cui/phono3py_script.py @@ -83,8 +83,6 @@ from phono3py.cui.triplets_info import show_num_triplets, write_grid_points from phono3py.file_IO import ( get_length_of_first_line, - parse_disp_fc2_yaml, - parse_disp_fc3_yaml, parse_FORCES_FC2, read_phonon_from_hdf5, write_fc2_to_hdf5, @@ -385,7 +383,6 @@ def create_FORCE_SETS_from_FORCES_FCx_then_exit( def create_FORCES_FC3_and_FORCES_FC2_then_exit( settings, - input_filename: Optional[str], cell_filename: Optional[str], log_level: Union[bool, int], ): @@ -397,24 +394,18 @@ def create_FORCES_FC3_and_FORCES_FC2_then_exit( # Create FORCES_FC3 # ##################### if settings.create_forces_fc3 or settings.create_forces_fc3_file: - if input_filename is None: - disp_fc3_filename = "disp_fc3.yaml" - else: - disp_fc3_filename = "disp_fc3." + input_filename + ".yaml" - disp_filename_candidates = ["phono3py_disp.yaml", disp_fc3_filename] + disp_filename_candidates = [ + "phono3py_disp.yaml", + ] if cell_filename is not None: disp_filename_candidates.insert(0, cell_filename) disp_filenames = files_exist(disp_filename_candidates, log_level, is_any=True) disp_filename = disp_filenames[0] - if disp_filename == disp_fc3_filename: - file_exists(disp_filename, log_level) - disp_dataset = parse_disp_fc3_yaml(filename=disp_filename) - else: - ph3py_yaml = Phono3pyYaml() - ph3py_yaml.read(disp_filename) - if ph3py_yaml.calculator is not None: - interface_mode = ph3py_yaml.calculator # overwrite - disp_dataset = ph3py_yaml.dataset + ph3py_yaml = Phono3pyYaml() + ph3py_yaml.read(disp_filename) + if ph3py_yaml.calculator is not None: + interface_mode = ph3py_yaml.calculator # overwrite + disp_dataset = ph3py_yaml.dataset if log_level: print("") @@ -488,27 +479,21 @@ def create_FORCES_FC3_and_FORCES_FC2_then_exit( # Create FORCES_FC2 # ##################### if settings.create_forces_fc2: - if input_filename is None: - disp_fc2_filename = "disp_fc2.yaml" - else: - disp_fc2_filename = "disp_fc2." + input_filename + ".yaml" - disp_filename_candidates = ["phono3py_disp.yaml", disp_fc2_filename] + disp_filename_candidates = [ + "phono3py_disp.yaml", + ] if cell_filename is not None: disp_filename_candidates.insert(0, cell_filename) disp_filenames = files_exist(disp_filename_candidates, log_level, is_any=True) disp_filename = disp_filenames[0] - if disp_filename == disp_fc2_filename: - file_exists(disp_filename, log_level) - disp_dataset = parse_disp_fc2_yaml(filename=disp_filename) - else: - # ph3py_yaml is not None, phono3py_disp.yaml is already read. - if ph3py_yaml is None: - ph3py_yaml = Phono3pyYaml() - ph3py_yaml.read(disp_filename) - if ph3py_yaml.calculator is not None: - interface_mode = ph3py_yaml.calculator # overwrite - disp_dataset = ph3py_yaml.phonon_dataset + # ph3py_yaml is not None, phono3py_disp.yaml is already read. + if ph3py_yaml is None: + ph3py_yaml = Phono3pyYaml() + ph3py_yaml.read(disp_filename) + if ph3py_yaml.calculator is not None: + interface_mode = ph3py_yaml.calculator # overwrite + disp_dataset = ph3py_yaml.phonon_dataset if log_level: print('Displacement dataset was read from "%s".' % disp_filename) @@ -1138,9 +1123,7 @@ def main(**argparse_control): #################################### # Create FORCES_FC3 and FORCES_FC2 # #################################### - create_FORCES_FC3_and_FORCES_FC2_then_exit( - settings, input_filename, cell_filename, log_level - ) + create_FORCES_FC3_and_FORCES_FC2_then_exit(settings, cell_filename, log_level) ########################################################### # Symmetry tolerance. Distance unit depends on interface. # diff --git a/phono3py/file_IO.py b/phono3py/file_IO.py index 20040997..a101c87e 100644 --- a/phono3py/file_IO.py +++ b/phono3py/file_IO.py @@ -51,28 +51,12 @@ from phono3py.version import __version__ -def write_cell_yaml(w, supercell): - """Write cell info. +def write_disp_fc3_yaml(dataset, supercell, filename="disp_fc3.yaml"): + """Write disp_fc3.yaml. - This is only used from write_disp_fc3_yaml and write_disp_fc2_yaml. - These methods are also deprecated. + This function should not be called from phono3py script from version 3. """ - warnings.warn("write_cell_yaml() is deprecated.", DeprecationWarning) - - w.write("lattice:\n") - for axis in supercell.get_cell(): - w.write("- [ %20.15f,%20.15f,%20.15f ]\n" % tuple(axis)) - symbols = supercell.get_chemical_symbols() - positions = supercell.get_scaled_positions() - w.write("atoms:\n") - for i, (s, v) in enumerate(zip(symbols, positions)): - w.write("- symbol: %-2s # %d\n" % (s, i + 1)) - w.write(" position: [ %18.14f,%18.14f,%18.14f ]\n" % tuple(v)) - - -def write_disp_fc3_yaml(dataset, supercell, filename="disp_fc3.yaml"): - """Write disp_fc3.yaml.""" warnings.warn("write_disp_fc3_yaml() is deprecated.", DeprecationWarning) w = open(filename, "w") @@ -158,7 +142,7 @@ def write_disp_fc3_yaml(dataset, supercell, filename="disp_fc3.yaml"): ids = ["%d" % disp2["id"] for disp2 in disp2_list] w.write(" displacement_ids: [ %s ]\n" % ", ".join(ids)) - write_cell_yaml(w, supercell) + _write_cell_yaml(w, supercell) w.close() @@ -166,7 +150,11 @@ def write_disp_fc3_yaml(dataset, supercell, filename="disp_fc3.yaml"): def write_disp_fc2_yaml(dataset, supercell, filename="disp_fc2.yaml"): - """Write disp_fc2.yaml.""" + """Write disp_fc2.yaml. + + This function should not be called from phono3py script from version 3. + + """ warnings.warn("write_disp_fc2_yaml() is deprecated.", DeprecationWarning) w = open(filename, "w") @@ -185,7 +173,7 @@ def write_disp_fc2_yaml(dataset, supercell, filename="disp_fc2.yaml"): ) if supercell is not None: - write_cell_yaml(w, supercell) + _write_cell_yaml(w, supercell) w.close() @@ -271,25 +259,6 @@ def write_FORCES_FC3(disp_dataset, forces_fc3=None, fp=None, filename="FORCES_FC w.close() -def write_fc3_dat(force_constants_third, filename="fc3.dat"): - """Write fc3.dat.""" - warnings.warn("write_fc3_dat() is deprecated.", DeprecationWarning) - - w = open(filename, "w") - for i in range(force_constants_third.shape[0]): - for j in range(force_constants_third.shape[1]): - for k in range(force_constants_third.shape[2]): - tensor3 = force_constants_third[i, j, k] - w.write( - " %d - %d - %d (%f)\n" - % (i + 1, j + 1, k + 1, np.abs(tensor3).sum()) - ) - for tensor2 in tensor3: - for vec in tensor2: - w.write("%20.14f %20.14f %20.14f\n" % tuple(vec)) - w.write("\n") - - def write_fc3_to_hdf5(fc3, filename="fc3.hdf5", p2s_map=None, compression="gzip"): """Write fc3 in fc3.hdf5. @@ -404,64 +373,6 @@ def read_fc2_from_hdf5(filename="fc2.hdf5", p2s_map=None): ) -def write_triplets( - triplets, weights, mesh, grid_address, grid_point=None, filename=None -): - """Write triplets in triplets.dat.""" - warnings.warn("write_triplets() is deprecated.", DeprecationWarning) - - triplets_filename = "triplets" - suffix = "-m%d%d%d" % tuple(mesh) - if grid_point is not None: - suffix += "-g%d" % grid_point - if filename is not None: - suffix += "." + filename - suffix += ".dat" - triplets_filename += suffix - w = open(triplets_filename, "w") - for weight, g3 in zip(weights, triplets): - w.write("%4d " % weight) - for q3 in grid_address[g3]: - w.write("%4d %4d %4d " % tuple(q3)) - w.write("\n") - w.close() - - -def write_grid_address(grid_address, mesh, filename=None): - """Write grid addresses in grid_address.dat.""" - warnings.warn("write_grid_address() is deprecated.", DeprecationWarning) - - grid_address_filename = "grid_address" - suffix = "-m%d%d%d" % tuple(mesh) - if filename is not None: - suffix += "." + filename - suffix += ".dat" - grid_address_filename += suffix - - w = open(grid_address_filename, "w") - w.write("# Grid addresses for %dx%dx%d mesh\n" % tuple(mesh)) - w.write( - "#%9s %8s %8s %8s %8s %8s %8s\n" - % ( - "index", - "a", - "b", - "c", - ("a%%%d" % mesh[0]), - ("b%%%d" % mesh[1]), - ("c%%%d" % mesh[2]), - ) - ) - for i, bz_q in enumerate(grid_address): - if i == np.prod(mesh): - w.write("#" + "-" * 78 + "\n") - q = bz_q % mesh - w.write("%10d %8d %8d %8d " % (i, bz_q[0], bz_q[1], bz_q[2])) - w.write("%8d %8d %8d\n" % tuple(q)) - - return grid_address_filename - - def write_grid_address_to_hdf5( grid_address, mesh, @@ -1566,7 +1477,7 @@ def write_ir_grid_points(bz_grid, grid_points, grid_weights, primitive_lattice): def parse_disp_fc2_yaml(filename="disp_fc2.yaml", return_cell=False): """Parse disp_fc2.yaml file. - This is obsolete at v2 and later versions. + This function should not be called from phono3py script from version 3. """ warnings.warn("parse_disp_fc2_yaml() is deprecated.", DeprecationWarning) @@ -1593,7 +1504,7 @@ def parse_disp_fc2_yaml(filename="disp_fc2.yaml", return_cell=False): def parse_disp_fc3_yaml(filename="disp_fc3.yaml", return_cell=False): """Parse disp_fc3.yaml file. - This is obsolete at v2 and later versions. + This function should not be called from phono3py script from version 3. """ warnings.warn("parse_disp_fc3_yaml() is deprecated.", DeprecationWarning) @@ -1800,3 +1711,23 @@ def _parse_force_lines(forcefile, num_atom): return None else: return np.array(forces) + + +def _write_cell_yaml(w, supercell): + """Write cell info. + + This is only used from write_disp_fc3_yaml and write_disp_fc2_yaml. + These methods are also deprecated. + + """ + warnings.warn("write_cell_yaml() is deprecated.", DeprecationWarning) + + w.write("lattice:\n") + for axis in supercell.get_cell(): + w.write("- [ %20.15f,%20.15f,%20.15f ]\n" % tuple(axis)) + symbols = supercell.get_chemical_symbols() + positions = supercell.get_scaled_positions() + w.write("atoms:\n") + for i, (s, v) in enumerate(zip(symbols, positions)): + w.write("- symbol: %-2s # %d\n" % (s, i + 1)) + w.write(" position: [ %18.14f,%18.14f,%18.14f ]\n" % tuple(v)) diff --git a/phono3py/other/isotope.py b/phono3py/other/isotope.py index f929e2db..275c6516 100644 --- a/phono3py/other/isotope.py +++ b/phono3py/other/isotope.py @@ -34,8 +34,9 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -import warnings -from typing import Dict, List, Optional, Tuple, Union +from __future__ import annotations + +from typing import Optional, Union import numpy as np from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix @@ -58,8 +59,8 @@ def get_mass_variances( primitive: Optional[PhonopyAtoms] = None, - symbols: Optional[Union[List[str], Tuple[str]]] = None, - isotope_data: Optional[Dict] = None, + symbols: Optional[Union[list[str], tuple[str]]] = None, + isotope_data: Optional[dict] = None, ): """Calculate mass variances.""" if primitive is not None: @@ -176,14 +177,6 @@ def sigma(self, sigma): else: self._sigma = float(sigma) - def set_sigma(self, sigma): - """Set smearing width.""" - warnings.warn( - "Isotope.set_sigma() is deprecated." "Use Isotope.sigma attribute.", - DeprecationWarning, - ) - self.sigma = sigma - @property def dynamical_matrix(self): """Return DynamicalMatrix* class instance.""" @@ -199,14 +192,6 @@ def gamma(self): """Return scattering strength.""" return self._gamma - def get_gamma(self): - """Return scattering strength.""" - warnings.warn( - "Isotope.get_gamma() is deprecated." "Use Isotope.gamma attribute.", - DeprecationWarning, - ) - return self.gamma - @property def bz_grid(self): """Return BZgrid class instance.""" @@ -217,15 +202,6 @@ def mass_variances(self): """Return mass variances.""" return self._mass_variances - def get_mass_variances(self): - """Return mass variances.""" - warnings.warn( - "Isotope.get_mass_variances() is deprecated." - "Use Isotope.mass_variances attribute.", - DeprecationWarning, - ) - return self.mass_variances - def get_phonons(self): """Return phonons on grid.""" return self._frequencies, self._eigenvectors, self._phonon_done diff --git a/phono3py/other/kaccum.py b/phono3py/other/kaccum.py index eec5ada4..520cf0a7 100644 --- a/phono3py/other/kaccum.py +++ b/phono3py/other/kaccum.py @@ -2,7 +2,6 @@ from __future__ import annotations -import warnings from collections.abc import Sequence from typing import Optional, Union @@ -134,37 +133,6 @@ def _get_bzgp2irgp_map(self, bz_grid, ir_grid_map): return bzgp2irgp_map -class KappaDOS(KappaDOSTHM): - """Deprecated class to calculate DOS like spectram with tetrahedron method.""" - - def __init__( - self, - mode_kappa: np.ndarray, - frequencies: np.ndarray, - bz_grid: BZGrid, - ir_grid_points: np.ndarray, - ir_grid_map: Optional[np.ndarray] = None, - frequency_points: Optional[np.ndarray] = None, - num_sampling_points: int = 100, - ): - """Init method.""" - warnings.warn( - "KappaDOS is deprecated." - "Use KappaDOSTHM instead with replacing ir_grid_points in BZ-grid " - "by ir_grid_points in GR-grid.", - DeprecationWarning, - ) - super().__init__( - mode_kappa, - frequencies, - bz_grid, - bz_grid.bzg2grg[ir_grid_points], - ir_grid_map=ir_grid_map, - frequency_points=frequency_points, - num_sampling_points=num_sampling_points, - ) - - class GammaDOSsmearing: """Class to calculate Gamma spectram by smearing method.""" @@ -245,7 +213,7 @@ def run_prop_dos( bz_grid: BZGrid, ): """Run DOS-like calculation.""" - kappa_dos = KappaDOS( + kappa_dos = KappaDOSTHM( mode_prop, frequencies, bz_grid, @@ -270,7 +238,7 @@ def run_mfp_dos( kdos = [] sampling_points = [] for i, _ in enumerate(mean_freepath): - kappa_dos = KappaDOS( + kappa_dos = KappaDOSTHM( mode_prop[i : i + 1, :, :], mean_freepath[i], bz_grid, diff --git a/phono3py/phonon3/dataset.py b/phono3py/phonon3/dataset.py index b6e11663..74fdea86 100644 --- a/phono3py/phonon3/dataset.py +++ b/phono3py/phonon3/dataset.py @@ -47,8 +47,7 @@ def get_displacements_and_forces_fc3(disp_dataset): Parameters ---------- disp_dataset : dict - Displacement dataset that may be obtained by - file_IO.parse_disp_fc3_yaml. + Displacement dataset. Returns ------- diff --git a/phono3py/phonon3/fc3.py b/phono3py/phonon3/fc3.py index 63e1e6c2..d7204971 100644 --- a/phono3py/phonon3/fc3.py +++ b/phono3py/phonon3/fc3.py @@ -692,7 +692,7 @@ def _get_fc3_least_atoms( first_atom_nums = [] for i in unique_first_atom_nums: if i != s2p_map[i]: - print("Something wrong in disp_fc3.yaml") + print("Something wrong in displacement dataset.") raise RuntimeError else: first_atom_nums.append(i) diff --git a/phono3py/version.py b/phono3py/version.py index 50d3a786..2a2c8289 100644 --- a/phono3py/version.py +++ b/phono3py/version.py @@ -34,4 +34,4 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -__version__ = "2.10.0" +__version__ = "3.0.0" diff --git a/test/other/test_kaccum.py b/test/other/test_kaccum.py index e62335b7..a88b40c2 100644 --- a/test/other/test_kaccum.py +++ b/test/other/test_kaccum.py @@ -5,15 +5,13 @@ from typing import Optional import numpy as np -import pytest from phono3py import Phono3py -from phono3py.other.kaccum import GammaDOSsmearing, KappaDOS, KappaDOSTHM, get_mfp +from phono3py.other.kaccum import GammaDOSsmearing, KappaDOSTHM, get_mfp from phono3py.phonon.grid import get_ir_grid_points -@pytest.mark.parametrize("use_legacy_class", [False, True]) -def test_kappados_si(si_pbesol: Phono3py, use_legacy_class: bool): +def test_kappados_si(si_pbesol: Phono3py): """Test KappaDOS class with Si. * 3x3 tensor vs frequency @@ -110,7 +108,6 @@ def test_kappados_si(si_pbesol: Phono3py, use_legacy_class: bool): ph3, tc.mode_kappa[0], freq_points=freq_points_in, - use_legacy_class=use_legacy_class, ) # for f, (jval, ival) in zip(freq_points, kdos): # print("[%.7f, %.7f, %.7f]," % (f, jval, ival)) @@ -122,7 +119,6 @@ def test_kappados_si(si_pbesol: Phono3py, use_legacy_class: bool): ph3, tc.gamma[0, :, :, :, None], freq_points=freq_points_in, - use_legacy_class=use_legacy_class, ) # for f, (jval, ival) in zip(freq_points, kdos): # print("[%.7f, %.7f, %.7f]," % (f, jval, ival)) @@ -131,9 +127,7 @@ def test_kappados_si(si_pbesol: Phono3py, use_legacy_class: bool): ) mfp_points_in = np.array(mfpdos_si).reshape(-1, 3)[:, 0] - mfp_points, mfpdos = _calculate_mfpdos( - ph3, mfp_points_in, use_legacy_class=use_legacy_class - ) + mfp_points, mfpdos = _calculate_mfpdos(ph3, mfp_points_in) # for f, (jval, ival) in zip(mfp_points, mfpdos): # print("[%.7f, %.7f, %.7f]," % (f, jval, ival)) np.testing.assert_allclose( @@ -141,8 +135,7 @@ def test_kappados_si(si_pbesol: Phono3py, use_legacy_class: bool): ) -@pytest.mark.parametrize("use_legacy_class", [False, True]) -def test_kappados_nacl(nacl_pbe: Phono3py, use_legacy_class: bool): +def test_kappados_nacl(nacl_pbe: Phono3py): """Test KappaDOS class with NaCl. * 3x3 tensor vs frequency @@ -248,7 +241,6 @@ def test_kappados_nacl(nacl_pbe: Phono3py, use_legacy_class: bool): ph3, tc.gamma[0, :, :, :, None], freq_points=freq_points_in, - use_legacy_class=use_legacy_class, ) for f, (jval, ival) in zip(freq_points, kdos): print("[%.7f, %.7f, %.7f]," % (f, jval, ival)) @@ -257,9 +249,7 @@ def test_kappados_nacl(nacl_pbe: Phono3py, use_legacy_class: bool): ) mfp_points_in = np.array(mfpdos_nacl).reshape(-1, 3)[:, 0] - mfp_points, mfpdos = _calculate_mfpdos( - ph3, mfp_points_in, use_legacy_class=use_legacy_class - ) + mfp_points, mfpdos = _calculate_mfpdos(ph3, mfp_points_in) # for f, (jval, ival) in zip(mfp_points, mfpdos): # print("[%.7f, %.7f, %.7f]," % (f, jval, ival)) np.testing.assert_allclose( @@ -305,79 +295,42 @@ def _calculate_kappados( ph3: Phono3py, mode_prop: np.ndarray, freq_points: Optional[np.ndarray] = None, - use_legacy_class: bool = False, ): tc = ph3.thermal_conductivity bz_grid = ph3.grid - frequencies, _, _ = ph3.get_phonon_data() - with pytest.deprecated_call(): - kappados = KappaDOS( - mode_prop, - frequencies, - bz_grid, - tc.grid_points, - frequency_points=freq_points, - ) - freq_points, kdos = kappados.get_kdos() - ir_grid_points, _, ir_grid_map = get_ir_grid_points(bz_grid) - if use_legacy_class: - with pytest.deprecated_call(): - kappados = KappaDOS( - mode_prop, - tc.frequencies, - bz_grid, - tc.grid_points, - ir_grid_map=ir_grid_map, - frequency_points=freq_points, - ) - else: - kappados = KappaDOSTHM( - mode_prop, - tc.frequencies, - bz_grid, - bz_grid.bzg2grg[tc.grid_points], - ir_grid_map=ir_grid_map, - frequency_points=freq_points, - ) + kappados = KappaDOSTHM( + mode_prop, + tc.frequencies, + bz_grid, + bz_grid.bzg2grg[tc.grid_points], + ir_grid_map=ir_grid_map, + frequency_points=freq_points, + ) ir_freq_points, ir_kdos = kappados.get_kdos() np.testing.assert_equal(bz_grid.bzg2grg[tc.grid_points], ir_grid_points) np.testing.assert_allclose(ir_freq_points, freq_points, rtol=0, atol=1e-5) - np.testing.assert_allclose(ir_kdos, kdos, rtol=0, atol=1e-5) - return freq_points, kdos[0, :, :, 0] + return freq_points, ir_kdos[0, :, :, 0] def _calculate_mfpdos( ph3: Phono3py, mfp_points=None, - use_legacy_class: bool = False, ): tc = ph3.thermal_conductivity bz_grid = ph3.grid mean_freepath = get_mfp(tc.gamma[0], tc.group_velocities) _, _, ir_grid_map = get_ir_grid_points(bz_grid) - if use_legacy_class: - with pytest.deprecated_call(): - mfpdos = KappaDOS( - tc.mode_kappa[0], - mean_freepath[0], - bz_grid, - tc.grid_points, - ir_grid_map=ir_grid_map, - frequency_points=mfp_points, - num_sampling_points=10, - ) - else: - mfpdos = KappaDOSTHM( - tc.mode_kappa[0], - mean_freepath[0], - bz_grid, - bz_grid.bzg2grg[tc.grid_points], - ir_grid_map=ir_grid_map, - frequency_points=mfp_points, - num_sampling_points=10, - ) + mfpdos = KappaDOSTHM( + tc.mode_kappa[0], + mean_freepath[0], + bz_grid, + bz_grid.bzg2grg[tc.grid_points], + ir_grid_map=ir_grid_map, + frequency_points=mfp_points, + num_sampling_points=10, + ) freq_points, kdos = mfpdos.get_kdos() return freq_points, kdos[0, :, :, 0] From 501f1502f303123cf1466214d9618d9e6eb353f1 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Fri, 19 Apr 2024 16:38:53 +0900 Subject: [PATCH 14/16] Refactor kaccum --- doc/Si-kaccum-MFP.png | Bin 24117 -> 24649 bytes doc/auxiliary-tools.md | 27 ++------ doc/changelog.md | 6 +- phono3py/cui/kaccum_script.py | 68 ++++++------------- phono3py/cui/phono3py_argparse.py | 2 +- phono3py/other/kaccum.py | 34 ++++++++-- test/other/test_kaccum.py | 107 +++++++++++++++++++++++++++++- 7 files changed, 165 insertions(+), 79 deletions(-) diff --git a/doc/Si-kaccum-MFP.png b/doc/Si-kaccum-MFP.png index f486f4ad1cdd4ef256f506d2a8977a07e85018bd..83eda0d29dc7f8ae53df8ec9797a8a591757004c 100644 GIT binary patch literal 24649 zcmZs@cRbbo|37}(DvInK$xcZ2h>Qp!E7_aO?CmlmqcX!WvyQz(c1UuJLpGJHjH8T` z{d*j}Kfk}eU00Xe<(${+yq?eJW85G2$Nh04)KnFSF40^b~`Mx*w*}t4^J!sJK`Dek!O~&gr6)2h+3HjZ@PM zF+Uj_ADlkRz5c@?>zwo&22wTFlNsJnoGxO)7Be?&TVtYwS$3_^C9EXnun=d! zi9Z^3lT)u8FSeiDj39+~GUGD7Y>oKME~W|55TWZ!SF&QA>B3!Re6G$juem=}C+iyP z>I$=;tk+l78hI;!87_HgT>FfzMz3VhnzxI+pX#H?ZHkwjrWH8V{!Q-QOX;RhD4K?q z@1H{;0?+YH*2HNP$_f#7P<`OceJ6o5=ck_H2zjO9U}tw0#~zt&lY-s9LqraDrI3xc zavxf*&a(-X$-=3kB$p>Ys*}kniF6OZtJnFf48*A0oatt9LPA2_J`M5Lia7Ak60sxE z`gKxGoAaFjgRf!{Demy0xw+Fh0!c!#47G@dd$7KO#Ei+>>D*X)q_b@Ff`Hf zuVS3gx$%XXoNv_NB#qJ8J*t{o^!gVt0@3camqKMkkvvJgyS|15CX=Xo&CeSw{lzn`KPGWdarZI!@CZX;}lrAFh0RfXW zyG@C9jw9w<*XIrK*o9V7CRuTYGEL5?3hknV2P}JCVFXNj-j8oyyEs{!8$azJ^!mY8 z%^HL62m8q|lL1_{303I|1C4ow?z`@*bHb_{JX=DS;FH^2MUid!v|-ZI!lFo3Y}q$2 zLvepJprQyjg!v*W%&*6qAxUSL1D>zlqX${A@o>Au88 zr=&#fI}t2G>TIB+H11GGZ5Ari#vAo5?rc~MDIzd1&!+ljmW`fc_XC+T0@gNd&S@@- zX!X(Hne^nXSjPsbze#o2{VcWPizgI!ztH{m?9v@Dc&w$R zf+<+DyUk&-L6XorKeJyozB+;bqsRO@b4-8K%^*oqYiE#=>+j@tc5F5lt(t!yDHZ&R z-GAczCD8!ZI7KD?^HOUAQu?=i zvvz_8qq&Yov$)1-2kZ(ac34a-sdo*FeVUliXPFzHdh(}q<$U3NBRd5zrK)i|!PIb` z^Jn)j!@bAthGg2+)1Pm3UtO=I3d=(Mf-`EIUZ1{HVW0X{WCw?HGw8RhisCbXrRr!U zeDw^HqBvahc59q&TYDpuMa->Ez@$H#8x_Jp{Z81`%*n`cjIorGNY6$vHJF}HP*9M4 zLZfh)mUofPjjgG?>I}juu=&IO5xHuiVv}M1ZCWj@q7lw%L6tZ^0j4roYoc6#4R_Iq z?Vx2i3RUfyX``2=Ms>hX8}*+dOM=p1?Tc*o;e$fV@_?A?5n1+1!=jfo-*P!~#7sBk zOSzQaT~uP#{K1Kt;{7>PA;iIVtD3R<&5vA6$mI&Piy~y1J-u1uvT)f((SF1B`UKKt zq#q`nJxsoS{kk-=MD8&4-nlqGmW`lN$Qd8%4NJ5Y#D6ygSn159+SuLDA7XsFDaI^z z3PBah$p!(vS=E{(T~NvUrI1AbN(pV|ZQ2fA@3`Q$|MG(Oj))}QixHEBSn^Zv)rxI# zbi^y&?BKN~cvh}(2OpZmZ@+$^n|L!sxcb%OQ8ztwmK;I<=6M|g>pJ!q^|)t~mE5Cl z%e_9X$Wf=n2z_QseP6zn0IbDgIAj~k?963mkGu%B3jvZeP=3*RrBNeESE#%>h()U zSAr*RR!)MeR;`w zFa6rfPVZ6Upvg3sSzK_(T?{9U8=ZpINTwvK@Spo4iP-%%BEj-*A?sCgWTw)^hpbs@ zWT)vys%Q);9m%DgPIl*uF~83K z9_U!&)A&47rM{X|6V3NdgQI^lC|`$@TvSmuh;%SPDb&_+KV9-WBhTNGLZyp^J-jE* zp3kmSsQH{ufk2%zp`e3`NFDs@Ihp8S+yfTEjt4B)ATq|(FKl|D)?Reap-@AV8d+)^ zep0XDr3j_mYa{P50349r?CrKSqTY3%D!DNN%fEknP6_l4&X#FP(;e(-d^R3zueIlF z*CNX`EE9`?J1^bHzF?hV`Gi6O^2Y%GdB5R8O(g06| zs6DD3;2*ho?9}GC=9RwE@C-*gH~K^?8SlJyG6cUaJX&FI+a>v05Ap=HMAvp(mTgh#s@5a16p62~=9xNlKxppf$Y7FxvO=_w=?alIIvQ(h z1@{hk7sag~KMctE<~Um95${}POs9eUdBgo!H{+Y@2#Yn#KzUI|J0c8j`t4_?OMNg$ zd-?vjg34F6mJAm|o+P&8X0{UTm4Uqds0wA)JHoi|hK501=3lauESd3A_6`mXj*ghW zL5e|SGQu%;?3_+5iX&Aqnt9$$s<)ikyc)@a9qbT@fQJ|XQQmADJ?l6%vIu=vVBYjCuui8piK8I}&G zxbvcJNOtl@4Pe+i*7FnJ2k(A##zKIv zy?(lBMw3N2KC8Bn`}Ps~&+>Lvj@tnMx$*J0XWIGhVcWa5b8S>{8m!kY?r}oK2aItz zZkMQ6UZ+ zWb<9ow8e{zTWyrbS&yEW$9J6LdzwgY&erB3vY1vGUPvlKmvS2&zvu-AlkzLZGkJgE z+4R-g1q!WgYDNDtdnb3Je-ok}PZI5P!RX#9e98nFd5>qo~V z%Nlp_{Mz|g0hOvj>$6f_bIzk~Hp~uAR%%SW-pygle|w4pUWXy!H!~Jn(m&)$ddPAuW{xFJktyT^=}I1(R#;C8()ZWei*0qX z9dw*)slnG6LjeJ&xym02-MUKlva`U9Dx5fb#5Ov8iU`45pZOk&8IM`NBvEsK^71(q zwYA!u323Ak?#Ij9dWmkof5;h(8hL&Xr8M&WU6`Ma(k`OBZ7A-yS8bD0qW$uy$&fQ4 zh=Q`QvGMkYjdqXbO{6}Lf}dEB1D7;?SG~)k_wlY|aZz_gdN#^-sbw^LdKMRCHrrBL zo44Y<;pV?KQIm>)JS7<~z?NwHqhpwYRlAJw1JYMg4pLhn)J71z_H~$-epLExV~Lib+c<0a%=W z8D{wajXh9Z()sd2O244FsLjz20%)66U-k#qYjFVzX7YUaL_umvb_NPHlr&f$RJP|% z(g!yM-mehly8kQu|KBglyHzG>>^0BY)&%)L89Fd=6|OO+kb$Z{chM_B$ztr&O%)9P z+bD}CY)-uRB=kiu!2L%s8Y&D%`=R|D@>es@8Ws{+l>u8ripJi2jeSIMc zWj|RJi>DY=BZF(=yqo-|M6aFTVEZ$8PD8n~Z6bpB$BX`~p}21f!(FiERx)~XwrKeq zseKH4XVG5b3|I+ELj0r1;oQtz3_N6mm~*S`CwU1u=J)pYR?Xl60^^phiz()6L}n8{ zR1;o|kj>WT;fmDYO!%{VKo<-P4J7?}U2IUGWZW){l|5#TDvT)v>ilVk8+Y*;dA#1y z%8FO|Hy(G8y_Dvv4>~odl zD&f9!HQO7;4AWPK#Tch0SCJ*L=ue(UYVYREPVd+7>)7i#3!LC|SfcWUJuqtHL%{D&rjLCMPqR zW3l^9QokDfi|d>L?G4%q`WOxH5_To$*~q*wt(8M+!!%8b=c`v4?LM3ei! ztTkgl&RnC?BiN0Y+Rf4NZ@uXJ3fQ5y=KNkC^XdKgI&g>bIo}%`jAY_e^U39G>L#mM z%qHChQ%Bqkw6waNg%{V{Q@`q+Fj;vvrkPg}ts3v!9r0W4iru*tZY*{_2G*^v<@Sq~ zy{#=nQigy@wnmTWT8vDV8c!0`!dYrg2BVe5#ni_#-Nl6`-CUC0&_BJ8FijA&XC~`_iO@6 znDWeHS5NmJ`_1-~DR~gJZJM}hG?$yYo$0!CITNx^CpCYR)6$}|JCArF`zq$_!yAuS zq+(1-k_$VkNIB+T*pBz8PZuqERj1y;;NWO;-mlNrw>zqn-kB(r7CpTwFo{Io7m0by zlSFN{?L95fXC|RTAkD?|pDuS;jsyZR2X7uc-BbF(gJYNAlESwaA;f<4imbEiFSx~C zYYa`z`ATN=&IkR`MSWS@?}X~H{v`@!U2zq`TMJgrdAy;tLeX*hqO$SEO0vsLov)Mr zYYVnz*A%|5c7hikDr;)l*XXrz63}mm?Acf$5L9`$Z+J29{=S{WBTr9uuLTx?GDbR- zk<8*m7J~~fR&W>TaxXj=X|SEFMWMQ|Bo z40ToP0u$w>b7(K%9SkcpJ~y9{1ox-;`oOr5s82P77`C(QKfiD9UEKFnZ=~i?H`rdQ z4t^$>dV<>Zr7gJ3xYI?0qs}I3db%@*K(^EPZ1PBK-B!ukxo0UqIWJe zO>AZ*x%NFcDB+m+=TVQF;d3sG9*U(#)0?xDqrys&)+u4bP;)x3zG!sx5-Pk;qk7zK zd%0hgOpdyu1VYbd1uWxxU8%=9Zx{L_gKW5fFjIR~oF*p_5`ME+OANp`#+M}=Tv&I$ z-8>{A91D(P6&Hj3&@Q45Q)ed&qc3d}p+N+4v^&!V!h^_6ifSlL*4}5l7p-VXMTuy= zI3Q+sHJ1hP-^S`YjY221mKOcj;j@18rx0DF@o$XB7+QJ1k6+AX1>olra~@%QQ=Oyd zX>=->VP^c=z^*xQskhS2$Zv7nVMr6=sK_Ke~V@|5FhE z(|Zwp|MMBx6W@xBi}7y|bfX1##7>X{3p)alvz>mb8Jz#l!g%z7nfiru?$YaV;cMUz}(Z10qr9IO2zDeuoIoZx(m>y#9;i zXUzgC!Ajm!h_EpdkJ&#C8XP^L(HQHQOXa8c9+6?}nL@DpCE7)|9}l`o=su2M(MKSL zJbo+$xr_dKfiSjhFEZ8A%5|{Yy0Iw!Nz-~&wL(|6$On;JwQ_{)>ryRu4ki^s6#aT& zsE@ zHZBdFjj|EjG#5SA<+Z0$*@wWeuqBujor2%8i z7OsN%*nY7X@w=_I#VxnWO>2Edo54HzD2}Hj!&_JE_~W<7&!zT?pa1&48&>}#h+3%e z&y1CZ)gUJkB246ys`a-FV&>$Bt6u5$ldjcN3FYOke`|V=zK*P;Da^Rt7e2TjMFj2m z4of42+q(Mx#bbDkR`o5NPK`%`w!QK;mli)j;Wa&0aP8c_WO?$>8GG!SquWeNJJqD; z>YLr={1JtIIr`y0+A4%x- z*O}BpLY*^Dwuec3V&^^*A)?kE!0p|*OOia(h2~x%V*WE~oLXu|Q$6fVm)GDA@MCk} zojZ@nGqrZVR&9E@VU(1nCa%YPfrP%5xLS8;sX0AT!r*~;?4J>@cKGZs6Mo|ti4Jn~ z;xa|JJCaRxS)OCv=O=IXOivb$VYU~G3pzxkt>?~zu^L~wQc2qFj3<SC6;=2T{9yjoeZEP>@q$(0 zsp~Qq#uc;FTKVX7@2#hn)y&kV?M_i%nut^HuKZl!wE}yqEkXVWDa~;aAqz?WwMI;x zA^#5PDfzik9-ZcGjrH#In)E~4;fmtCndsDqJqiOK$PpADrQP4C=4)x)BDccR1pcl^A zXRV?*@Q}4&QKD7zQ`6)QSHSzP_jOwb`{ zp+p!K!8P{);=JX%gSDgvv8P&eRFakS&bV7x)|9%_t+w~wPY=HGsjIr(^lHDtYEx>@ zZn%QGp}EipVm*TE>i9+3*4c(%)h&j8J4MO2#AS?%1k&JR6S8`CfpdT8a?c~l2$?v;YUfSuxdlj2R;E`ChZyTKkujp|h(PvGD>fMdId!M! zRCzBhsa%R!bQTOB=<3@>BVpHhZZdk6dUHUSLr`>f%V0}JY!}hUUdyWBlCiG)x2*Tq zxs!T-d-~ej-@=#p7@r82^am_Jk+VZ}q1hf-@b{?>tVIo;mQDQbI7C;c`QDuxNAog$ z9M-o>d|!WiDCr^P+``*}`)3fk*S|kwtI$(tFa;qN z(aN+FUNpLrXmM^|>?E!A_4ScFT7C?uIt#QvP*6k}<#{PL{XXogs*(L>r z+6il;0$9}a$WbuW?>529!abeYFUJ>d|6*s}Xxg!FaXoxg(&EK&jDEVEsvdp@!7Bx2 zXgGG27$Z0F_y;m4^BSel~%*1=?}0+moZE8&AQw5LErI(W}XMieS+VF>_96 zk;Ry~pU>-OIuoR&-$GEl`qm zT9=Me&z9}G=M|bbsCv-_tFi>rii+ibpwk9Tk_O_go55i%{UiqR@-2=n_1QpHkMB}2c{z~)`!FCdog z{&iEHXN>z`>FGXXTjM*cqctMX2dwQ6?g|tNGA*zyT~v!{pH;>bLOB{rF+O-2|C+(A z1f_02%H%srDoXta?U&xMHubXLH$1s&W0mEt^4AqNQSLhvo<5=trBD_ZjRcXsgBA!% zsW1W`M9xQ0A+L%m@e(rytk}*!y?x8a_W5*6)%aNbKODz@5php5%)JuQamTgO9lWjv z8BRL6YV2cB@pIx!@}#|#fm~VU*Q)t0^>l7#wBm7umpf-0Y+HtOhluD!&-*sHg!8$Z z%xyjTiI#dNP7KK#!KMdGTyx(?BQCs)cU{&U@A)>PCYfOsrNxqsjxonK4D}OGRz|EdcT*Yfb|qYm-+%|8$J%b>3Kog zd6tWdjnpl;L~?(G;Bj2$NAhq71VVdc;!R?+#?@oXDr&duY!9lVe-∋CJ11N<8=- zQ~)L_NAwb|m##xux0Tt`F=7>9lU}IkNI8nnggOj}Zy{8o=U(Z(4~t<9BmHC((6ptj zS|3Bc3sn+MQN?D+qdtq`8oSr{pKO3OV+vZ-kV>Inu zug9wfN|Jvx<=d2vftXk^%fbzPONE}Ecay`DL}KRDzI5nbr}lZL3oQ~e>FsbuM|bn| z@!YLnpY}R2u`S#hG_W?2(oUFxD%8* zSgATVoo9XR=PDk4wQAwr*2zn_)-?u{39ZPx8wiBcw^yqy#!o1=rBBF-W!G28VsTvG zD3uvTMJ4FnPdqYJIfDpP3%6p^Q0m_V#YTC4&Wbnx+calmVm|-1=PStG5vvx_1f59| z53vLs^_y#;tq95*?=jZMMsbdG{UGT7yjfQ&*iiN;Jq73<$M=8WqydGGHJ?HR-heVl ze-tP&SCG<(>J-~d(47hX1W8^%-s%nIMbQR_pA3_qqy!UM5pWCy_?n5w>`P^_CTHm0 z8`4mk_81BI@V)$|k@pcczEGmpW*w?NucvEBO}jUqWXR5v&UpK-S4V>5T?hZd4dK1~ zQ-GDBIMYl_Q-+4v(8CK}0-}?mh3FIo;`fi@(*e^`HaJt{&&2M^G zEZEXq&|Z#tL(k02&JdX=A))^Oki;u~vmL|c zMg>v(!vR}CCNcYhJCtQ6iS733HJ&aEL!}#=^|zx6LZ;zl6aB@2PXW>d@j}hB zOm`M{SEYAxGdStbu5UiZ*}$XLxvje1vAOp5*v-QqaUBd!(=GQKMRfV$gRB9bMODlY zAtjegT3oH)pkF+XQ{6p8t{bHF!OL>Q%~d%zG>g8v&KLW}-tc*RAUZ;@s|6_NcYzqBLYHop>e&+tnwncLN#a}un~-%hLVW$e^HTH; z__O8N)K0GXN;0O6tK+eocN2&r7RecuKl<2TZqM3O^j@w$A~Hdjo%9 z-_*W~iN(wnM~67ry^Ml-%zN>Vv2Eq@nky$Zpx%`P?;qTi8#Qsd+!TadR~~RIAH<+g zz)XLA-R5~vQWU+u5*+T@yzCpxx1+zUK_(598r^+Qche`@OK(j$a_vC_LA1bAPH+M= zuU%AE;*wAIJh48gg_e-X|Dtd8XTc?hiB!c|qw1vQG;-S0yBUb$|Ub6jBk%XB_# zvP#PhJaBD%+gDbAfK;}|rx3;>p=FS

-g0+Zj!57QD8`#OYU}0yQ2$YF~(5?)b|U z!DxU3seC1kLYbn`Kb&TLOZMBerG3)(XTo6>+wgLad8LrpZ}EPdpO%)J zLB7Q`rn8Is1X(9JzH4V3-kcxq z&^*)l+J~?kzvA50aNxG1_v0syCZf{Bkovt+MnZCb@)5n zBP1E_@j-`s4Mq6WL0t>ztqPrI@8Y+Z+a!P{L#DZKSU~5JCl`u=c^1by{#6~CVd5Q# zrBM|pY%qvyTa_W)0L7P0?496J%KmxGQ8ARe)*eHW^}K#{U($jlKGZtp>%MW9ZqW$M z{+f0$gjm{GgS`us*Nh&vG`;p&<#=|`&j|a}29{B%8Odi!C2!QN%Oop`Tu8%FRK#Q% z;@gq2-J@T>+D(3So37;U8jGfdR1o)0lB>sW^KREzf565Tlx|iJGh|QfpmxK-F43Y! z@7AarQ3ruu>RBRn0hOGIk>B_zIv`MjSEd052j&9AfZ9df-t-(6Mt|P1D`6XYZSfYl zjh>%v#%&@e_KkQ9WEpkAq1c1ZQN@j#5w4LxEzy2l7(XZ6jJ;qmQDB}>8Cfwqgy;UN z#@|t9HKh&`X?`6Y?e=fXU1&P=26WNlylW*SCT3j=Ud8sYKy?Sm#60SsvKYO$#4@lu zsk3T0GN!ZB-_jm`ZkMEqpgS0Qs2)GCuX8?na5wyzr$!P z{?_|iCID%Z%h%#FskA{xs}9L;Bc?`AMQdbavP>kt#cdmBIdAehtfH%Hyr3Z0*50NrJ@ z+eDdOboe*vQY^ly78u9syB6Y-u9`otutFVt^+l(tSxau7O`j#kna=yNmKIR6YOZhk z7pcvF&DZ5d;UhGgIAgWz&zZFV&?WrdZozFHi^tHwz6$>bjJjS)7gdg-W0j|l3g*}L zUU?Bnhn~liK3JbVN|+4*2L5=e|G+3CMw9}YZq{IM0>2zItb_y9|bVf<^tY?GZ4*Kz(aR%$Qw_H}#*bIIiaG&(lKvCy`T zs?g=m!}r@aUp#YzW@EscqZO&Yaz`(Xye1B|F!(wjQSKCG$)@qyRa8RPHG9RoMo%t= zKOHis)?!n@=A)xk%40MJa5S--FnsJwYioB?sxa(6K#G-{G6Ybzi@$zQJ!#Lv-t9Ks zvI$9g_8?Bm=LqY|fL~xM@BK$)`^}C|v>}8xYVsu6PrB1tZCr~2Sg|Q>s1OoARQ;XK|MZ`cJew?NH+VO%Qm@cB|WnTSFyY&`=)Tm^-^X$ ze`?>A(%0{yS;6Do$y}fNre@@;9QqC{N*~8J<6+o$Sc<&dQ1k|qzAuS_PLH*d@rAz z+}Qv83J6zBDawd}k+RukEFI`KB{x~wKs zzgZYz*mx5^pk;8x+${w0jmMsUeU}nr{~z6vV)_`$3J$DFP+@E znsBWv4`tgpy(139hB{GkHap~WY}9W>`NE#EfWH{2gWW3iZaBVjQTgpii>SAG5f*RM z;64HSriV$#KcyZX{mwt}edAn0iK;Uf{-hf7pGMr!ox+*`x4nm< z-8ok-K<>NCBqoBoes}kyQSqHFbiq}5dVFfNm2rQr8J}XknBkej8o=x!`S0dmrL~Ho zuyX#GUx@O~Zb%nC8oeB`ShM-a&RYyReByW(29FT9R<4T`j@WjZ{(kWY>UyZ#SV?X% z{B9F8Q*QaoVzye+)whUB%8r4mZvZvlbZ1+31}J`JTi>V{d~8ovKTIl-~X+;8V{ zR0Y(K9UP`$*&Tt$IYk;dgr9?WT291)ra`5q13@E z|L%ICMIUFj+>rNT<@Tx$wON@c*00``?_?KLLpG4~dh|fRjoWE3!=wiYcWd;t@SGY@ zd!xKL-*sqYS>+@z|wxojk!Vfsq(=rvR^9%4P?qGKvrNcEk-NV2!1Xe|p*4 z+Iwj4Mm+uqo99B>4CLdmK(qrU+Biu;ORKjMRDHfbUMgYRdvbr?I0+L@L0}v=atR%{ zBB6v1X4KM}PwN*b{xov*D|Og^Ejq>bvkYD!6SxbB!PNlgC#8FPP@r$~!?$0b^I4&l zU?qu@+fnaSJV9e=+~3C`v#68ZeEq<@>Iq&;a7~6P?7MJ`300a4MqrV4QU7kkY&HX% zt}yI&a?1uvS4(SUy8$FVl)M@JqqcX^A|)y8k(KTD-9KKK{Porazpc1p?@J19u%rb$ zYa!6sPp-Ldm^Dtj8`Mp%mEy&lx$Inx*=oCAj}_?0fJSK1{VV(Ld4%s*#YuAoge3E^ z;vDXDLD0s&XoTjo@pmX8jF03kAkUbrL4y%sGAoJPXP^D2&?6$ejF&SFQl}utY7Kaj z!VN%MO`4caz{Htw$wCU|WHS0MB>95Ems5s7OY6?tFXC)$LJei38CBM(D-@t>Lok!1 z<0}cCBsT-bL~3BB6lb&uKsHP!2QU}LXxxD$8eIX&WxWa2)>qd?h!=}5Ui~!|+wEba zHbZ{ihf;rG1IS(kFS|J21tPNnVN+|hfR4#V4^SE$6g*V5?wdYK!1oVBWtfnvleH2V zs~>&#-8M0gu%CSh=kOMteWzNdIL`*iSG3tk7h#Agh=`0LN!01hjxhKA`z!OVkMa>Z_ctxb+oZbB($^tC^@UlqUENs!j>0 zOFyJV)VjIPOQbOagcYV>lARS2j&6G^>*=mhs7Yl83Vw=16SE4c#(;g+_e~WQZB(0J z>*ZNW@dIbRpirvMrmnWQKWX{#|CPD|=CTqnL(= z*f)gJGh<_!4?!gMe%M-oOz!P=-lN4sF2G{YY|mr>1bF@V_wfi6O5-ziQ!WaXXERZ& zqIl8G3?tAdpdzOZ#6`V+-oxs);YoRScj!(Ogad@rJ$vF?zx@K-4dhWHC=}{s%%*eB zSxy}_!fC)sb_pnVm~2qSqzPeB`Ah{Y$U7BKgH@F_5FQxOg<^Ch+ic$)kIG^%`4C+3_J;pdvwbD z#^W?E^u0tMsN&plgUIYD7-}Q6nP8{#$4RNuoiq<%s_tQ(( zBmg_g)gbfBWMzU1=V|Y4`g#wal}F9TZ->VV@G8B)GvD2yeuG2rb|&ojH2++H{r#wT z|L~U>C?H4!Dj05O_u+xxX=?B6Wb~>N`o2%4%Dk-ePrXr|m5c-ZVRbN~RzKCPh#z$= zK9~c#T70OX=#7w6_XSpdt4R1E}=e_9pnH;dta_9(1VE~@Z>H`{_&kwBj6j8aIbfrX$g9Ba#n!SU^0s1CURhT zt^rbLo()>Ke7s?3kND)KFKr+&_g3nU&Znj3DRxBCcb||CDeZSgdGpuyXDum@e8})q zq1@R(4i^+Ef|xt0*Y)al?WkK5j;x6P;N-JKHy_J(MuAL=nwC&LINHGO{9pIYW{JduQ!`?m(#1mxf5S!)5A+c)ES z=IKd2A5S%H6tRoK#_OuxuOj_rF0~4(X8obnkFyOZWB~KYHX9d!f=QnTV}n#?C7qZs z_bM_yEGc4-zh{Z`73m#xxcMDbI+RJKb|KjsbSsqGb3O+H*X zSp#^CXn6@0nCang$G-#5AzU-R9YU~IUP9=<#^1Gq---9Z4KQGc^^UN3Qfu!aiqz!1 zt4mdc`o-^X=vrHPG||~Y>qZ7|^%I$ckAjLvA-_1{vwwRuz#eJ9B^<@2?-?b`$)hj6>@1FHnG*(!ji4g_&MhN4ZQQQCn=nCe`y@*9$4Ul}?!V0bz3rjrp6e~A`P zR&(Uyp#!9%b|+BM4le99g>qI)NrlR+#Dm970FRwC2r|?^yelzQP*csiNDIGsF1+&= zX#vy}*Hx3hkSO36AgYHY<*OF`+gyd`%$0Fp#gYcYLHBC>B<;njIj-v}XuMvW9Fg#w zk$zR+(c(XX<#4~oi-3mQrtJiuQNhS?6kq>lVrh6&iyOIB1ngVlgTHtJ96!UjHO~xW z7k<_cL|j&kpr0ua3GSkw@;9yx7s65YyB_P?M!1XxT!WiIWiSbB$uQ6h&p5?6+bO@H}sxna->OX}r= z@IoxJHFG=ni)B^X=737Zxe!&n|39<>i?yvu>9d^s{XWkoK8Gr-O~fk7CPjX^SMvdf z-z?rE8TvOt7{^VqQqm()>M`3G`aR&}QtN53Ry?!kHjB%R*HA*&XS0Kvc-&64V(TL> zVFzyzV<>TiEKEG1>+L5&y0D7r4e+-jox_LJMD&1z-M{(8{2HRJRh06Tb@bNiRq^qFJuN+@>6)T59JA+Z6gs@1Y zYWQ>pa$s=Yq0VP&6jxX3zpxOhm;+7_SWl}({WF z=qPc7ymy}K8gqP7s0pysKGmNvr4Yo+y4>=(vrX2^zq-2Qtsg*1HP-lPdQxiZZ;sgo zJOQM@_t))JBD*N4m}l(7X9vNUi7d`l+9bkrx`jVK_`1N>&=?gol}-h7txg z!)*0b#MN3x4TfZ1P5ddhY^;97_p=|=sbIdC^EVkNS_1P%6DCc+T)ix0^x~kR|KYoO}5TBX6z8ihVD=Vw7ywP}$4mU733$;A? z(f4yi&rZB{xtSD3m4>co81MAqwFJefnab^uzl?mZMVsOS*YN)&1Lv-(U>3d-Sr9_f z50?@r|CdAZ_TAlqwL8N9k*ukhWz3#-cqD0W<=(((C$ahSfB8R(jQ)Pkf}R)j%4~q6bxrT*Bk-hRlE#%8qe4yb{-4D3$|RKdt+$=&~*5relqyd{pWccUusMT-*uTLHowFnsI{RVE)HOwdL%d?d}n zBcM)&{~ZgNVFrGD0;yN3PVr-3(jYr%Tq&2oj3Blfuv!l`Bl=(7o>ebWfT2LAr=#1U zkeh+2`IIk7@^7=6fjq(ppuL=o6-M1)khDnV#Xyz#%X{n7|MFe(18DUrk*{~6M% zu@xJ|OZui1(3$qJq(w-)E1s|=hJd8#@(QU7Nm|?Wweo)eoi{TYi2v=!>jXU;kj*hR z`m_GSO3dG_-1SF6u683mvB#XkDNAmBy7^sz`2W2+u2?roFQI6~w~xd_nJJHJ|8gej z#o2KVvP=7Cpsfp4KfZ|oAy4Fjx9lU93j9P-`BaR$2k;s;ke(hZLFBTe#Ejvo?qMd_ zrq*{c2ILP#E~#A57nnU3Y{ZM~@Nqu*ldzRLncQm)Rf+h*R)TwY!xK^kiLMTr#jh;o zD<@S2*Y$%hR~A2v*oGSZ0c$%vC$9R1YJUYgW9~gHjlVgbd4~+mO%X1kFuJW#C{Xl* zb>r_`tUNv5-Ic*yA-ur9TW_HDLZNxHFQ#=K;Ni_Y95DCuOZ-u=X>qegvHzZF_}7>5 z*FJT8r(ti6&k^>iiuc1P+OzvG#@f#g4*TsVNORbbs$7x$cIm3z6saa^X7K)m$N%$1 z{0C;qZ?C!gM4%ptKGg|$?8ozOe>H&1*HAJTO=7>xWKD*zk4{4ktfdu7ycvy$V}TVW zcoqilOk%UftbcYq9RIflX%IA&u7K=3wNH?SNCk7A@s%=bIQbX(c}#{B?LT)v2T-Kk z>4#VPqF3-$c8ZrT6HmU)#TH5fy6)E+dWrCj5;BbCf#<5Q%^qC8V08%NODCfUFnY8| z*T*91PJRI&o3*7+1}DHn8{huxtt@sGrPsUG z7zo8k(u^eIKdXZWiI@{EUD<2SP|ddzIyZpxL)T;XgOfstyA3|fPe!v0ATQVa*fwtji!+O?JaiJZDjFuQAJzSwh+zb_*Ly0!+?`lxX z?fg)}*EZYhbm}PpAZg;Eu+q;C-)#=CKCZhgtd1`>uYoN()b)xm{;pP>7^8&mP{|tb zxoRPkk_Q5`X3BgHDZm#3eI#WZ9*$D)GHLNZnUmoElq6*=$aEK!)Z` z@PthR2I3i87Kjy8pi3*utX`3dSzu`xj&PU72K^%AeR}xY0qt=qpX8 zgzL^I6TqLQX{Zcg3jTxZFvlHuKE2dyhSLQ zkKg_p-L~$JrNMMB2O$-R?Xa7=XZ|)B;-CdA{)Gg$~B_|NqKyMfmy) zikBVfpsXs?q}jg(67wZJf|Bl~h8c%Cu7pdA`dV5_Pp{Q~opfDNVYhin*C9u@#q!cNRIa>Yb zHti{%-e-S3uZ~!3;r{{u5Qr0j869lNA`|PG>YjYLWd=PWqlsZ>h)2IxX=GE_}Y396>>f zmzTFVZ)>Hw*-=HSIDM6_3;7X=n;D27TWqoa5`C`dW{79vnUrKOL0=R(*-Z#wq-cvr zSs)#XDm9-BhO|G!YifZ>s$@_Z?szC@^YJ(hqINjv>7yBLLLjM6@tfE z5MExy%@5?m`IVZWtgmc1(^Q|qCY@%!e{|P^WGD%i<1|S}OA1_%;g{k8TR)KWr5VEZ z+~l6o@-1EPG<=yt3s{{L?%I34fc!n6&y!+s9Zq%&)dSy0)}TT+^mi@)XSAJtOO@*u zU>%U)$@s~=mc9QjCft(#dlGZA<>q7QjA_rSwG5g^-`d6eA-3SZW7mZJ9&~%c<=x~# z<3jB3i5H5C?;h`@PsHrU?a)&(M;OyrfF~K(szmZ(u(b0PU&wuv=&pl=;84aUnGQ5( zCw3iV(@aPef4J4s;?;Y5@vh4s2{8lHKbVZR|IZ!my?NV%3xK(KFzg%8CB`;b;F^yC zxjBCUJPK3DO-lOy7C`UpNbM(=xih{u!!>8GUfgT^!$BV)@_+H%Z#BI%lek^AZNz<^ z!SLC$@69{8Cik|w2w5n32IUpKzez7HkMx#b&sfE%`~?_!^b9c4wLi~eEqD(w0>8Gp zG&2{o--SuN1b!on1>{wHVsL1Y)|`D1NvDw{4l4d^P~USRqtWVI5!iCzI|K1dOPg$1 zl>vDbOlqlgf~+7uOA+WPPOQrPP@!B@xVYGxdaTOpD916e3_|+!bILL;av7EmaR*36 zlrQd+CcSvPKFM&17oU~t*2hl>{QbF!XRu6vNKvc6ysFhtS}YUU8@TjpsTStP3Hm(`YS`<&@>#z>a#=Ofd%6@+M_V%yvnhb_=|xDj@DS=Q;nubg<6o z>~cit#`gBf_)lTHiGacP@y^rcI>vm4gXRmG`nnJ+CuT7iF-y;uWsp{dm`rA@ z$^*jnAMCYbw?K4<=t;WI&?;p?2)>=pgYE6`kPy-$4WnXS5G&cBq}Zb1B$NPre7z+97eyj+j*t{sI%+Z| zJiM)=V~G~2I$3qiN>f!;)zapdvs2jQ*6Ql&WJkovA`Q7Ib4xH6i&cDZ(#Apzy%M55 z-Bg8~4^icAk!#vw0?v#bVUzPnwY^7bWYR0aczZUwzgYI}YtdDQFdKd1><1|4h`Jj2 zQ`0jLISP5&4`MFy5O8AcNUo54?F7D@qqLGNS)P9KJy&^CB;j5q5&8h&4dg9ZDjFLp z5&N&LkwOSRk?_Y-n>&sdtn(TkT$YxK$F$!pEE-__4?Pg{OpR6KM5=ym2IwR7@h#6; zAUAUn0i%O+HJipHDkMae9Z)OrXq4?+{^Tnf(Y%#X>n4PA<~jH>F(C}{eC90`zrsO! z?xiDj*(ADxIt}U?)>O^qeXojg+~r;y2B$n|Nt>v;21+U44|rV{V%2G{uhV+Z95Vi? z)lZ1E^_i*d%F0SS1k+UAe!jxAK1ys5c6Q7Hhf1uD(yeX*jj*0}^pf+GhAFyCg!o}8 zpfXpW46@#HJ8i{Y9aGiPohPaGOjufD6%mJ#!nc|XK{vkoz)WxDP3k5B!ttY9X4BI5 z!JtvMPD_@ud^BWL2fbqJX?x**TA=@^tP+frJ$VZWzsoC%UTM*q9$F|Jv@3>JjNxFT z_QGkG=^Jsu+(h4Vrq8D--@QVdaVf>0=&eeLL3Cr5?L zC%F+@ zBcOC#gp;+;lbpdIYDvaHnAw!F+Mki*0Frx>eBsmCgVdVD#6&8}$CHtCK0^!OGc1{t za*1%7Xs&t~u2LeBVCeLNrX6s;xYb^HM4h)Esm|~uZ;?J)-byt2NQGvmXT?DYN&?&Q^qXt%}7=Q(d}sLdss50l#rWV|{U< zB+3L01IyR5sXe$k+LhqOnW+CjN7G&ot)pk@|UY}Dd-+|P_z5q2(9AAo6G^qv`1&F&U6fN&%8f$E&XS$!`L62 z%6ji|&rx8y+`RtBF8hyp_W!!%NOR(W{LJ5BjqC38a~$hAS28Eocw+}3f^N}^&DpS5 z$Hb@y3$I&&7DVFwqZBzlBc(+)Hhpuz$E)|;+=U?(Az*!aPBjg%q`zq1;=5- z03MIeOwBvWtlKp?8RmetU(;4Cu$S7eKk{tx6@rf|5x-em9QMo7@oh=JH5W1 z0b-dq)I~A~-uVRxF*PlzNLK$Rxm0>IpP`YeX(ugGIX-`I4j{42X$cC1hljJrPE8+w z4@<8aIpyC}7^#qghH?hd8$WOTfW=ByV`6T>M0*481_xXA4JYfcHy=aehO%Zu6Q=;0dxu@!i?%XMpvLaeVtgU(PSVFTrY^MMMLhaOLG>zX74Hc~v& znyifzOqW|pKKz1m>DLE0e@+yA*+Su=zLehf+-05avw6Nr$i%CIRFsKR_ECmU;!!{n zZ>Zo9saq}W-!2=zQJt2}_GV6|0Ld3G_Ow@;9B2$6MqjKLjFWj0y*@2R0TPOLa{%8dB)kJuoMCX})-<;4#FD|e zo_l)tvfIyc>M#FXS>eubWwLGnP*#pCOV{89N)J-31*jK`2AV9sD&wQ>ZRA%LjI;6p z+hYM6b`{q}Sw3ec!;hV5KO>Ts<=HE~wo@=TrvT>APBL0k2mxkHsvz~u-%P>@fni}b zqM_Qc8(16JNf+w@l$y*3X0ZFB3R++{$N|B3S3$zeBQ#byI}PXa&T27TPM&Np>4I36 zNH2o34}GsZT`t3?fA*n|e5;>4B}g5Ixo^pj*g3|XV^;NWh;N2@!M5+QWQ)TtHYjt# zGFr0YZYxJ0OcH-GiKz+qS)JSSs>Rnd;{^8?l$z2SqtU38m`r^cjZ1GwOPZhL0tgol zR}59b#iW;wM@}Debkgds760hHrXEw-yiXCjSw`hdKtwPbQHpnaXaQ~5#{^}sh$AJ4v-UF}fstL>N#}ui?RohRs;S0^W zA~ud?pB1W~yG~+}507|JoO0yt{_!ojR}8=(ZxGit3J}MW@On=x){CD9RQ(DPfPvB!Ms)5I4^pwk4dteur`my7o?0=1g z1k2~;o?c~TLqhA>LJ+r;55@6=xGUNiULfL>nm08|FLg3o5wO1013ER#3D1|M@=L%zVQg0zGUp@dpjn^r9&9el3xtDz%-EI(>J7DUKk<8T!_A$0MX;(1;3US zaH##5IflEmKKi#Kiys*JlO%iWD63Gjo%Wap<%C?(8Z4bKz24Jf*i3%=SC~=BnOO4q z^yzB&*`h0aOS9N)slR%Tm2u*M`l-ChstXn3EFAdO`INxY&~F-OoNWFDcvJt%T!He8 zFkmVwa`JyQu6kdnyE>CH-3ICn1AQr>FDO0nrIPGJyAp{2S-PwC%7{U)YoQlrghfSF z)dh9cY&R3>?BY^XnXL2oj2SxaO^za|xU7_LqtW*64W28F&Y)g^7At96Ki0XL==L3! zzWAzWAYq6%)|tO>mo1 zc*oZFlWXcks;6CB0s!$4`6mLF;O1e9&ecS0x=W^@y_)Cg(V_&^zT=SBs9>TcBdpQq z5&NYH%6VbXl;2fCCYHP+EBb3;VL?_a;v?|B;AEcRdu#pN!<(Z$#i+cHw<8UITRF@| zhzp=0dfGZVDP|d{_3-qgAW*TsU(F(Az-PE{I`L-VOJWWeYux2BM8}UP%tu=qzh(SG zm2w|cso_?3oDwq->_9f?qVYl- zdHtjW2fb##@XZ|^h7j=F3USSk^$9=@%|gX9tb_RthX)4I7Sdh^xc?}1Wg6TBmBi-! zx33U3wI6xA4u?rGd!5B9{Gg+(-$=6GoqVCBxD}Lf)DRIQtfoHr>G(#LhoU%hcjqCE z)E?^fkIhEJ;&*ko@rPPU;QheM-F{fn2dR$k%9z~%u*zHXzG%TEXheXV5=8#X;0-+V zt>Kf;Wm8NEq)(%9woEq4zc$n$&R+86?HGA3-uaW7WkOI5+m4H;G{VlFTOVyyElDuJ zEDy`LR&tJ;s)x7fQ~U?{h9 z)6T|e_$vLmpM&zP%KpkufbK$Cr2btj!&7^$vTwq9;BK$d=g-CR-h8W}?Xy7z?|2$8 zCI39@&E4fJk9~%85LIl6v2O4%CLu$7P6zwnZEElLPf+l#dUkBSvN6g&E-8tRv%rk2 z!pez?LA_CO{IP>I0&tq6)u~@^=kJozT|TsRpf4v!mb|c4f|?y7Tv@p)c3!K_QdiBY z`{%NZ^}j~c>+9>AZ8(4oD+zqR_dT@zW!R{x|L6HTq6AO{1t885e8LI$Y-Ba?x@yv( z6cc_54;9TI*xol=<`{bko!H_Jgy=oDuRWP!HJ*tgBC|H*gTvCfkcc<_h>rO?3$XV- zEl7$CQ$E+}FR=sy0gJ^t*XOmzw2&ily$$(R#=(=+rMjk~h!PDUPJqH1=q-ep1p+g? zShF<4fBh!Uq1x4;mr7xB+%xIBAkgb3^l>)zIoa18H^ne>I6?RdJpv%hbB=D_frBAN zxsDvtH(U6?!^^D8XXum9oB5tGnh7|UbRnYd_Cvns&Tw9b3c0G`yqYzIIF*PP(?N zZ+R}!J0u$y)JC8DtoU~b6}JDaC;))YD8o(P1~MWJ<$}&+5P&#=gE9hb|HozbB77)B VC%#B+7sQhPWnyS literal 24117 zcma&O1z1)4*DbseMFa$tZV{zRx)GG_?vj!Q=|(9538e%iC8WDk6cA9lq@^3_hC4au z|GxKqzVCkb`Sv-FoU`4mwf1k#-<)HPG3E+Ul$W@RMvR6)Anr;@K37H{ZulV(NPeg{ z;GHsw=d17o<%NvIbHvs4KWU9Q(Fnu?gw%5pRku$YQ|`K|i=A-)*t0>C?12*`QLBB#QFL8H1eNnJdbP_yB8Kr zxVftryW_ElIbL>tcpSut1*bH`z-^hRwxba8I)06Oo7^&omY$1iWMDw~MELrj_t21P zt|~DO&NEa8XJ;}_OI1=UnJ)7ch3d1Q_yBNJ0)b+tz0{{B9EcfI7- z;-6QmpTjj3S=6Q$_up{4I6a7nh#EyN;ZjWv=lWLJDz@|=L+wK%#VbG((66$u4}0_vZ=6STk=9T8-X zOF!wDnX_r7PrmN%?zVpW7Or-^5SN5Ze!|DI{#ZiCPfxb9Wa2~H=Ioj;swyjEGrk-j zAJ^2>eEReWuAup7+Os_m8-aN0D{J+bQAcM2hep0YtGvCt+iKy*`>V^NW~zq|-|6AP zLAMg-m^^b;6O)qGx3_!E$~j3#@>96?^RQ$*elNELZ+`U{Ufxurh+CQ4OWZjcp78k% z!NFJ2`(Uu{)?GxSuC#R5)>OTWjLh|w<>lq6<*O6bKP+_8M>K!=vhp*5y|S`$oV372 z&-<>q&yPRff+Hz}r;|_lFZt>mz0CD&kkV9DRH8q9@>(6pIy^k2;cs^T{4-9A01IKc zIaWHK6|djb)wR37-{80$0f#DO80casZn%$!m*tVGiy=%)`&K@uSi5p91j}dTheXN= zWf+8L?W?RH-py1Kl4(96WzTkznB#1AbnV~^R=lBE7&bZ3BQ`tHn&(wf4;&5+3L z287SEXWy1neSH?YV)@L_wN0W0D+$3 z&*eo`2t~Crx1G6Y8inkvEHw>{iQm6-s}@!TVM7{jZr)Xby_SQ_ub6$az{hMwCd`kcmto`YE-Gc=B1i>j}Noulk1$bm@@mP#Umk%dKrEAN8(TWf1xetwQZ>gC1JjO)R=;u)vV zkgCv0 zq6xAtdrNY|Oe@3r8fRx`GtIsR?y_!f$C1L9UleKS>FGB%HUt6^znb+kDlfttW&qwVZ` z@Vngn4kg@DEWF^`tYeRG_EG)y3=jZ2UKR@>GR7?ri%)^p6dV2CF{kq&J z?@dfhjJd^NC*pcZ9otb~TYGW3UW|3Y$CskGIePqeTUjp`SM2&d1Fs#ij<}H#bz$LV zn~ADiLI2JHA^&$(*tZc+%^@woD)KGRg3yUce}))woH~~YE%AM*ce3%ZvBOZ0Uni}UgFnwpq+?e|MqSXfX}QU+iP3k!?BjgD54mJVIh+u7NHM5;0I<=od@O&lK& zFBenMCSd&iuKav!AUXlF?r~ov2WhdovvZ{`^VP|6YAlOHB>BaNjvdjGAwNIAP-?LU zlpJ1{4b93K6K;!R;(96Oc!d;!-|g)|z02?k;jIZz6pYpzh(WQO!CKv|507O|4>!-? z>aay2<&lz-qF|CrqVhV-1-v?oTOG{VJ30UU{kiEp?D^HDrR=Ky(b0JM917l70obqh zSCoeBA0cC-1u6(0sEUiDP(q%k%F4=8$(HS3|9TsaBqKu=Ae+R+#Lay)R%#3pe0jFr zMK$ab-G#$O_-<**>~5-9Qt=Z`&V#>yzCm_`q7SiA{`_I?mw#TlxvQ5_9w^tDGPk^> zWUSV4`RzTbCCKnL;}uLgRc@Z16JulVNO^a|xu-R2Uzuuak3*@xXhmn?;p2n!b3Su* zL6Izvg^s>HDVZ~PeOE{IJ#8R&J`2X&6cE?&_44vkFVsGm3uDpN(z}KzBa+aaYy)ZU@e7_UQ z_|V~x7o{o3RE}cF?p~xf z&Mqq(n3_r!|9mH~b#2X(f+96SBF_dF^`>&J>ZT%{ceLcATeog~{rYw4x14(Q?NQkN zYHIL~v)5$)VbY@j*owc$$LUk;UEAIoe5PyB?px9pKp;MH7C^Z~7{>i?*o8AN@A_;l zZ>JQ~r-0olXevtf=5?e>`uJQzRj0;A^*R0%9269EeC%edkr{iP{Ez57Oy-Tr&NaMXqq{3+Hk$c(rj;RFQoAX)+h0E5BTgnZ8L zvZ2UXq*cCjyfY7x|95}?QGnOqueT|D&v(LDnvQ3``uq9$K|JvB^LKpzzA;trdbsi1 zyysH_hs6NI2N8z&vUq5QYN-h7d&xe@r7b~Y!fPrNQ6s*R0}bL5iLwY5FgM+$%Z z_~Co8ZjhW)C$xJk6-|9lR9ae^jh$U>TX|wkt^Uwl@g79KTz}5UASNc}yW&Ho)bkFXqN98C=zDAHi~W`UysYB!@$vij z@5eFejMTn*JnlQQz9zvZ(Rdp8^;8rAC}8z#Jtfim+U+2dABr5002qHPN~s|Ea4+G^ zR#Ffz1;T$-c`Il?_3C^%*>yTMKc6aq_UNyav@|^(9WT`8@bJ5o-#a^VOG*?EiU5_J z9c|~^c+E{kb#-6AdQ$`^lIP{*mZzd52_Q49T5@MXoC&2iL-E!B*mM)ncXk2IW{`;+ z+@bG(@Hr44>rzTdKK!}9HRkPSM274re)`>BuFu#WZUklOz^8P>ZX*!o5J7%J*OEGGG5srj|(O9m?m_sz1bVDFc)h z@4wDy6dzHhw`gSfyl|R$RU1H@nZuq$&7LA1<_>*C0JcSvs45M?$d3SLHl<`;RT`Y( z?zS6&&mY92F3eHGPXct@57f0aQnqcSCa;(n2{{-d9w&=^E@vV<#E-*B#OoONAoW>C zhjoLIaU}p-)Lcii>}3EpakyDqr6GU+rWawAJ_fqDoo;n4JL7kwqBdlNQ?#yZrAFsC zedeZG*$CQd>QzzP9jk>_?SwD#@~B3|^6M;Ri%1p)#IU?_^aQm6E{RT94?BxTfyan1IEb@)b%|bczy`^9Ixw$6?qXw}o zuAAcrhZA;KFeu5v3>I7NmEl(TU#5n2qBS?+htDI0ZOUlud7IF z!TbRgtxYcQ!_q@3$B^f$EY&-PV;<}qr)<`Zgd?_aNkKLgL3$J@6Sw@GW#2d5s3|RV z_!e~6^Kiq$(h_#gHdH)Q(?!U-P>vv3juvRmZ)}`vPFc@1d0$;#K%sp9{(Tq1*3NF( zo<~Yb3V@zyx;<1HcC(+FOuMJ}wl{Mi(}$#YRybCJMEBk=Z?=nzUtW8+0Q3HO-w;j>Lfd^(#u+~teEdexkSm_QrE4Ti#9j+sPZ6g@U0q$4 zy4PV4Zw^?~dI7v=E8nZ!DgW=aVP2Oki?d?Sd&HBz?Qoo_dj>d{tQ&z_N29S5|u-w0;)dS zZl0-kwd=Tqw6)xmeAYqXTcB63uHju>QSrU4Et$`GEm}cHtK#LoaPfJ2p4iyfJmi}< zL=+XHg4iE8AU)lt9dOdn_161XpHrIk{Qhy`vlV4>mH^ZBU#u>shCIph;%AcTI4z~6 zr8v~GnhkC{5oG)^On{-h=fk;mb#T|*4-YA>uDDR`puZPEGLez`MC`7G!hHWZu*koI*r z#WCxh&Rku(0!joFA846ZbO_6qlAP?l*o7DJTvF!WD<9P@$I_3RlL#0Th5-Id)D|^M{y!YLkcV^ zrhhLl&I$_(isl&YZ>!&3~--6!0 zMV&fCM?6J+kGfEqkg>MDsWuug%sq~$-1Yd)Ch1p9ui-zx?F8&Si7=iHa%BJf`ST7A z^)-Jom@OAnAVq(D_|vm9SMrdldrL^jPmiP~a4We3Ik{uwYcr~(=~x05`95uOwT?U2 zRO%T#Viq88Wc>a=00sjH>$-j<(DYYd#H)|T9(lZX}nJc(R{X54?N4}+Fig)NtaHju=RbP-4l7A zNnhE>SJ_CfHujNW3ATHC;htiZBhSPDS_nFyv02iS6?`aMiD8e2E#OWr+MB+QZ487Z zt$uGkdT&JDh>?$^9@nqwf#0u6H-sD8=5X2gRk&J_G@iD0zitrSHhpu`ef3b}OBq7Al~O};la>w#ZY3`sqH zzxzY^>(?(YA%i3#1fou1JmbqR+>vZ_s{ygf)z*<$xhjO+uw=d%%d5qc8)sTX^-mE~ zEt|v$RyK?L46b>*kaeDqi}%0o65HW4b#+p{4Qb0slknKbb22hA0-gjKhr23@M^MuJ z5|kx#Pf2%DdmV}^xd~)QZCTIPVX znoC7BmC2GmYV64P+^&Fag6>^g$=92dtact4vsG1Gn8A(!7=QebwBhTcxXJ3bx+3sO z`h;|p_tPeB*DDM&tup&?biq(^WKBkVejc8fxVXBSnwNvGF9AKilT@WC{&j!EJZ#N2 z#L^B;XG>~_l)wY;V*1VH-)~qzr~n%yVAz6C`wT^t>tYrA24Ram?SSsX-1^6P%@^i# zp|{YMDH-)qyXj&UNx!;ah|h^`XUw>_uq7};G6h=An2HT6jT@?NyYBURu|GfHeHi!o zu^np|K`Zg!`^&k#Q#?NEH(#Tt1u`S+(&FyAH6-)9$|rG!%U?0%Y;J5|s7+AKq7gm~ z2B_ftJL5XvR@E@;4Cy{9-R`hQGiKPP-xJl8W#-DdeH=z^C(>=g2E>U_>h{7yrrhmQ ze*{XOa@)4eW11w$&ueAR@i7v2ec0?+yB%1&pL`5y^ha>?ZL`n*O5%nx=?RSB*47rs zcBzrV$@O-EqJmqa zio3e=(JSeGMm!u*;rCZu(M9V-b;j-D8vNq!t;^{?dJ)j=&p| zTvd!za&U=WaCvfYncl-O!bXfV9FN*tO!q^2V7=$Qhb$~Ci1?%Yg1!sFDJ3Ok-03qG z>MVJ7StDIwU=)~dQ4DX4OejMqlwkKpn zA2%>+Ej^u(iB)25A|E*nI>8m+QJmw7mn8E~=n|tyrW$Jyku9jeNOL^(Mq!NaX3?v2 zg4089?CN@yz-}&db-qVNu*k}M(;mydx|VrS;T2N^Z7%~uvYz3fB%R6qOK{pRPa zmRAhhxLjUSOH->3nTA-}qFLEKbU=Fdj0MSSeT0cQ9+;1rh;VMZX#t<}V@T4+UQPc* z@R%$e?V6>W#a!bjmBn-Q1h+TASF{AFPi%cIHszdJljx%E;5k=6uUz%jUqVBBEUAl_ zACN(?=u}yE$1wp=figrW>~juXxSS(B>y5FE%#pZF7Fl6^mx-xyAG1`E=tksk?kOAs zQK}L{4~n6_n4zk5zu>va1}UzqZ!xfhRf?dcX8(Q*@Zit`0+KBA`9*4B1OA+lyfH5@$hRO_#r zWT3-BOFbewANTWp=d>EOuJ<@Wu|JARealACS4Y)X#^h_}MXiP{Y8XLR*->SN(KlTL zPB!fDfT+2-IXEl?4E5v3ZBtW%^z_n6J!$b#^a~O)<=iv!?`antUVrv+TQVi7i!BXf zyCZip_w$$YFn7Qk(U50;cq{>JBtLf(%$9LmU(_2?()&$|$@u|3`6TRn34JGYpzeQX zzq0oV-~gpXMn;B{49R<>yQ#po0CSr-Ni4W9-nH0ls(nba-4slpT#6uy_|?$pts zb`ZgJ5fQD<&Ao^@F&!zOFlafuHq9UyW z>WjMPC3-<+dJ-jM{vSDHba!v;lb%jR{Dr#RMbZzn6Yv{=&EfUm)q^DYO@2AP&y!aFg8ZMQKx5QJYR@bz$D{) z6yGgJvxnW%?p&C?k-Sc`ijPl1l`Swo(mS?D_(u2nuxS`1|NOxKvSNqWLB}0Uf;5f( zdD<1wb{17<;T{7g-T{3~8-W3-uulV5vPiFeapCF|iEz$v>=liGM}&bJ7P;AD+Yh5N z_u*h`b|IDnfpRDJgAokzMkBUVzbPFeU?24bB_gi2>bHeX7NSE#(Lk&*PxV%>y(wXD z&nh>?b|fCaVE7HE)6`($yy@}L?m$b91}1rT*@ne@$oKHwaoP8;CYk^Qc*g#ui-FL2 zg=Wlq*N8u`ti+cBLA$z$zdzmRNg?3I3cCkzn3lG-rK3xuD1v0{{}4E_yvFm`P7{d;($20fAb|(xUDN?4XlwI1tthI6X)?;o$sO$N?d|Rk zOn%9Bx><6m6fhp^d-3_5fdLt-wVq}sXR^^oD0yUsno=rNp31aZ4>s45lyxH!fU9NH zDw&&5f`7eYoAMMu7lq6tvv4Gq`bom>K& zDd+~In9C0QpcLxM$^8PFKF=%R(k*s$$yaENa!(Q9Ioc zoetb<@u14VHCWO;XIq{-rF?IPNoryAVt$)k*(EP)VL9$}BIA_wUU&Zo?vHKgM*o)W zErh`LTZl&1KK9weLl*`H2EL2V;09_M8sG<@YkHGGYl=wabrer1!^g)Dhys2KT0q$G z&5rXRef_fnkpd#WtTcRH=&&6MT__vq+%}FMd@L<5c>MX|bevep-?~pHqCpWK@J(kM zlPQArX5)5HriqPBp&Ug()c@ebVwr$-L_$L9N_!k|e0FBg=%HrwIcQLbgB{Ji;PEXn z@ArFW1iS>JSRuN&tOG;T$6e!i2fQy9Nz=y}m8cA?N8_KEWT~_WFf%(sYqN-!b&4TF`({bvP_6e*Es$eDkYPt^6=nfolXYe z;%9!bfaH1KkR8_{+hK1PdI$QjGNxnt#Cm* z|MSehgBi=GT~onJ>|j^CsfJPmuxef}{^$&Tf1yPKcAN@r{fz9H(bcOzm=$Gd^U?$e z4^%m%Pf<~{VPSNk%9r<`i6I~--`d=qs#T(B^p_6*;{0d$bL{4? zcxo2RSwh=hhS~tj8YLN+zbkZ9iBSd67L<)R*9Y2X`F_M$7+;|$9wMtYL z94rA8jWBeJz-&ekv5OzPJPaw+DyL>=SK&p?udE#Y_HBMpK{z}tth=+5h|_WibOE63 z$7-KEd?@O@N*It;i7d9tFzDm%+){GzETe#NpT)_x=D)lEA6;TT0VwBK2qoG3e%JN+ zL3P_NvHPCF3^~~NECHR@5m#a@@XPwc6U4sRDOrH?_H<)ROblf7ru_kcj80ZlrJji~A9c|=~q=M*kT=6C^QKfH$wheTD`4oZ9 z7o{L3H8eJYW=~8&Q0=s;Y{NG~e&hJg58P1B)f`{O&Pn#QEXiJ)<^-i(N10!8Mk*;a zPZStwV>IgruQ0|DaPUD&zmSs}L|^$Tw;zx|AETn6r857O$O$>^x^SEj9z)5LiD#CR zkr~W3-lX)S7@b99vaO}%DXbAt*oS>jS6W(H1U(O54l)oi7X2Su zl(dWNocWtxFFuzzy@7Zsb%AFOD%JCc?(2oslarIG;Q;tSl={SGiWYBmzx$EM(Pm>( z)vLab4~p{RmM!vQrY05n_cDs*^KKwkehuBW2ZH6QM)Z^Et(o6rXss?_w+uOlBM^-d+2%i>sRa=ws9>qo z5L?~#TwiT}cl$Qib6F3Y7Gd(EMy#SL1E1yBH6mMwF_jajRqLo$If~G(&KCmx%#NlF z3>q#PSJD*F@h3ShE-wGHO0yE7%LY?X-j}oNPmHI|ghxoS3PZ!f=owfl!_OyeT4YFZ zCh+bcT8^PU0hpFiP-yl#VT@i01@- zW=?i2FP3TK<#@7f_%J&HF$7%2=-|HwE3Y##$YMKeM zb^iOCm&+-CtNHy0m8=|XmNvDZhF>IQ15&g%?SQ<>%UW+wsGZ$Fwyt#2LVf00iq|9sg~Bu=8xUV@c;CjTV%HpoKRQfC+g z`@Vf3*n#~cy3-_@S`G~p6ZAJnHMLl-hh4tPQ8`z*~4orIchVWvLO8WzuYk#7A5VP z=n065i^CsQkTa*mRIF8U?M*P9H~rRcFKbo1zgDlNmb*zM#Y`lVyY)eYIiQAz(-_8rnqj1sM zo(Z>!l2QlGsFFFmjE?v10CLDZ^b3utl5vPVyeJU2bEw`vA`E*4_8JrqqFb8_`^)zw zT%a%x6gTxol1fmQS81m_Es|i7m6r&(FJj>3eGBnwslo z4}zig(CwO6!m~u|W}gJym50YBCrL?2bhNZA9`fHm++~<`jJ@T^LOy=BFQd|cOmX^N zRC?+%oAP3CTqKN?N}XS)>~8}B;%UczC9bKxz5TOi&p@%xUpLUHVJCbCL`uDC%OkO&)LYA|qCkeYlN}*Yb-2S%PT*l`YY0DWLni;%O98bjv}{?&Vgs ziR{ZS@=MCfSnl1P*DS>pAo*5hLQ6uj^|rXctxuiK)1|(?qnxX~5F2sbvM{ZHTz4A< z#YGo640~>E&Ew+K8M-x_x}n%4Gn``x^#~46cU=Qds)$bH1&f)xeHp+*k zz3}pDiyB0w$m$M$MR0I1cn+=&EtX)V0SmK{wgZ%VY-&u4Fwzy3u-X`F4B+7n~ zys9{~|Em{N*!HObiBiU0(g1;Q%%p39ltcs|dT|jC7gtqN^T+aD^X2I}ICv(a4lMde ziu@F2mYY3Ja*p5k_U_)H3_uNiB&C!kZ&dI;wcg>JwGH%zm-|9!bqx*Y|9S!B<>f)) z0sVb2n|*3+2kG3jr;VvLKJuNnX@H}Fs3euaxueNyi8HeZI@votV#HIsoEgM)opZh} zv+oxNsQuT0BFN*{_UiYXRfpY3Y--|h-K)sYU12un`pmwT(GxAT)TJ0ih^H*@iz$50 zdCwnG21sTN0EY}+!^+AE_yTQC2h&rCTc9ag?)BX^`gO9s5^uet{zkNuVq0T48$C7r zXode)H*C0%HydZ(Y^|*HfA&1u5<{%3nR&Mfs~k4Y5tRB!X4Z%}>&>g& zk1qzN$Hbs|RIe-&l{|lFVP!=X5ECCy^Z0RKKmZ1XkN`V7I~CRZV!Xjfc6o%QrRCcC zI&^a~sDF0cz|m7w#LUU5MjzY${X4d(uAZJmIB{%LRG%6msq*&buRUgpIF+HMQ1=5G z)845=qaDWZYvUqeL?Qx83m6T`?@)q8Z4$b_2KT*~&FjSe<`(+GE9en1GcpuXg+ODv z=7Kn_M!*3C?aeh<`-UZa^X5%W>B(%I4-Pu%T6McMttd%ImnSb987?*|Q9~x*YcZ zPg(UyV&*mma}D_0u3bTu0%sfLqUnkn6Zg7I*qGxz;nsi(^XAPB(YH?@&|aHHuhzK= za&pW&lnj&mFXd?hiCQi*_BC!}c&eN`bv1OYP)%1&V;F0IQ#G-VumAFs9tvyHucaNg zYp+G+qeqWaH8j@28IbCGdENJ3TVk)y9o^k|qdU)+lI(a$wG}xBs)u+|{E-(9~uFD$wk`FY*<_IoV2mGd^dVuv!49`7DQMg;G-tAK*b&dx54C%`-}7J2BEJKr!o=;I zuR-s=S^KBXa-WOSM$eeW@`S0{~f zWWJlg%BJ2XTJFDUe!{}ah~0=lEXg%niEwjs14rBooI;o2`xqxqEX|ONFB==y{QUKB zf?VP(WeZ!|qUL7dq(TcRG)(8e#%Y)y8^?SC~*z!6R z&$a5{jQp&V{8#FetPwZU!W*ud^Zyvof7k?EULxJmwyiq zn3$Pi_l!H`GGt@7CPr?aM$ska;UB zaMSsatevV!D1J7dQwgS)GI<|oaX49^L#Z0O4(c{M;a^EFnZFvxw zP{Ea8auWXiPUlQSBW?w?yY$Z^5qzibEC_@^5OcAg1;;oriO|^H0@p4OKg)}Y;LkkU z>lFq+ZBRk)kYd1T?!w4gQjWNTH}}=iV-1Okog4ust?LuJ2}ITe>M1Gh%ts2x#l?Y9 z5oory%(yI$bN1OitgWIZ_eAj@&1bSj2~xM3mROw%*|!9c5&NR2WVht59SGo@oSvQr ztD~!|W?! z7MGXf!1Ee z*%=M}pypeVnxH2{V~3d19!!S)wX(DChCcr7+qbT+H83lvR)SW+#l=-ucNiMhyk2<} z#N**4YalJc%;NcX$c#G^?$yiUZV@#g>Y$1M)YPiB)l*bdG%})0?Y#nf?{DiW%60a* ze9p%TwZ~4@f3B+g-^i2veW;!YS)0u5JD!uf`^kGMiI~BtY=u;Ee0;N#E^?{j>a>nW z+nV>ze>!@qkzdbMb z#UuE=fA3bR;$k-_jL_%~vK$jv)5CHYVYm1W(%!pN{ODu8B)!`dLlx;!W6I1ozj2m=meN z>{lZAJe>|u_Q)y@BwZ8ps&mHL`6I=mZTnw?4c>c_PeV>Z5hkdRJH`TgpolhBoNN&7%?-3^9Yh z{W-6l1edP3zJAK?46azAZtd>UuOxDE@^L4NGTWgk9RjZR*aPa67shJCj^^G`O#q*0 z5>bqRt9k)#fGqVO^bzgt*)=uCU`JiaU1i9@O*p{2C%_{59vz)tG=9g8vOuC4dhgi~ zCtj?N8DEfZV?vE16ZCNU{`Q`jm>AgOeP7xIr7XxBSH|Csw`%K|)@QL$2NuL+BlP-g z>ui1`#KgF0f%WyCii&-OI@PF@fHP#}__F-{y;pi9VQs&7${wzF z4_4JaF)m{gyhN}L=xAthj9P5gVJa8+f;_n4(#AZsnF;Bv$^dlw0| zqRFZ^cP32#0HUeYEwUsM^r!&kGpgvuabOHZE1AF<%RW)V_a}x;{)aR-p#&+g?D?sI zr>4$kf(vEd4*iMrqx&|e?3@kOldIQia{Cod+kIt3$i6h2H>UvJK5>F ze!8M=y-Q(ZqxbROKZl3D%%`;f12%IsXwYc1sPt@xR2)_(|JiSK7C!S)*G_>u=W$0+ zA8ByaiGd-%q=eO`r|360PgC8Oh`_W9%9D&tcSQvUXu=>&gHhDLUITbdZ|@c0DDW?3=aciFHL-4*R{!XSJ9BN>&eto(M>^ zr;z@AzbI*(Q7kMk-@bVh$XwcGPEKxanMWUBzzA4o&~tP&HL*n@o!Hvgz;YNJWLdeM zu&B_Ul3w&)EjG;l?(#mpnoKDJhr5z8Cp0=~G#$#)%7E^6Dxbei-Mc zh4k=y6zh@iZO$1~!bH!+Cw?E83Y-EZw9>X;ie(qC(c;ACJ)Eb;pk6@B%#06XHCO#= z9Z=u`qQDFsN%-Ui42pqzP; zUvqv7;U|%H4V)8~kbJ?!|Dx0lZVMJ(x_5Cv6u=PH1Xa;$=!>F9Ht9vH;f$C4PI0jg-f$wS^+Y(Ro-cMC*Tea*fN_>*BHPIrmpTX z+wyuj#nblx1XoKxFN8L{wjh69{9+M;k~M<)wnIqU9Nb+1s3n@Gb>a>7e?(wH zJSkfSEA)!&FX;e8{yWJ*+Ly2a!A8fVSuzXW)_BQB;4Je8U(i?jK^z|y=Xyy#+l=B; zbSQZGaw>>o@Q#6?c<ZXh7U*2pMU!LdzxHfg0LPY+g66_8UhOy zY25^-A@M>k1?5FkU*88hG4R`gPta;Qd(<#xs&34C2nkN1hxv|xni}qMmuwkf8h3Dkd^lIeQ;4D>M&qFZvsZedWElaAHK&H zr48kC`)dd^8q8lIbq|(09-(=K>7vND-h2f2o=xs^8sIKr)278U@M=k`wShh%1DnF2!`5Cv(RD36P*8)r!J`L09%(ginwl6{QC8*U&CyKYaD zL4$F5Xv~^9Z{aG6-r&ByAJ7Z)sY^xO2wQ9GYqJ&6TDB6ET z*z`R4_;KV+@sBoN4C)7PKcy5M(fOE}w?IKsYZQQ0>qt#pJ;dl3a1>`aXG|Eo9yXPI zlNeh8Lt8MHQSado7O@nWwW2K+plZnTyFyk0PQ@DhAh z(9@IoI)P>a?3|1YyUvd}@nJ!i=;6WQ!zV73aHYP=!f1RjnDp@dwYFhvYn#^&wn&iZ zVx!C4Ji^ljVR$=-54F%d?)j$Nav3uH{G zyt=*D5*4NS@WJD-Q$&;FGA89CMn*;khVGJs6@3eaVlYc*y z9emq9ibtAoJu*A*Az8+6?%644(^P{hMDns3*A$Zbe|A)-)bJco zEK5r;6_X|B8x-@imqK?d1Snbr&X6g>v33$VZ12mH<*wx>Trn~Mx03X9N>$sO51-ZZIm2#Q=ax8E z4)Ih_4oMegRyM|-B*=0BD+eqYs7OiO*ee4gDCw%FLz;HR1n|K${p#!|Qdt843>$!$ z0wRYk;U~h>IXH+{(_pH`%ZuV_cO;NQ#<~^0L|x8>fT3Kk-yKS{z&RihVW;!^ob!k& zLCx&TE-x}MHxFx@17`}I0gV2_>@y4=(&+C{J=yLl|4wn;hf&K6GOb)w#DjzEe0+6H zO+H|O1w1@aVG-TK@TI5-MB+J^*NtU~j*f=32VYX?WnbE3nK*&3ce2pj6ZiF1DSVUT z?!cAHBo_IBpM=|bjMZv{(bLmY&{;uA)!MosYcbV0_SOvmKfbR#mzuHrU-f0iKagY z$jr``&Kdx)Ci*r%fBlQnXK?wGbNiK(kk{ZjCNGydYtleQw9T!(bv?Z%AufJj6qwZw z{(hsG|J!^w$EVkau3ha(Vn=6O{AHw1kYPf`$$uo|a9P#y*7=`ciqN(fqLOB2*}JW( z6yp;Uz^Lv`xz6mk38)#v&_l<&qy2qV^#Qy}sy;Rf)pyG`*Gp8P1TFG)gdY+WQ2VLz zKL8IR9T!*Cd?ZX`~o&m-=z@woV@~J|<{@BwqFsOB)-V`O< z`U8FB!1Z$h-lV0abqG49WSYe)WFCU9h**g{714xV{(i_pVfFMYM@JP-E^cn%0Ai=S0oZmDHLXlC}QXh7gT2%thr0Z|#j7v_ z(yebxKBXofN~NJ$&=+N?PU9toew{MFda6#1IbL{w5J+~b!BqXJcA*F(u`HjcuVGfw5$&%7(=M z>&n=dPajR*bhU-|ABz05_-Q(Iq_?IO)n za8Kyl#{jJM_@9@vHOH^1DB|y%1uhH`1cU$7eJ7dt>NE)KkwP7qnIDB+GWnV^Y`rOz z__ZEMn9&&_x+km%p%N5#43ux}?Pj|RotT(rthH6=KI{^fHXFb)A?joX9~8hdAAons zsHxGU#SOjeh-z~-FtrIo26HvAL5Qg=3t=N@AspzW9}l{tOggYHC>*9X>#O^$eLE2^ zggxvL7D2Rc!g05?>(?*N)}S~Aj&KlXvsI;u0%i1PyK#Z;4oogjHn^+Q(xvBGWLDQt zPyu3Xv3^;113#=SPpJ=TG_=;>l60!et|knWO@H}0Os}TuAGDTM&yhtJ{B4y7>mxj# zLbvuq_u^V}UWW>UWhGwem~Ql0H!MsVL!u}H=zZ?DWR+q90|SeTS;E@bRE;9jrN>^( zw1QWP_`22c+fu&5;Nj=@Yz?>r=}qLH7d&cZ4CF4jbmwVrIXSskuXdq(7uPp7Vc$?_DBHM(3?zL$LuOuJ!x7e=ImFv$|5qSC%K z!QkOA4BZ2rAO{0p;GORJE+;P!1N%BnUMFcH{y;hnB=4zKSC4a+6>K=%f(ONzJ(EM& z0pUMZVg#|ce|V@~1#U#&yp|9^!HGk4z-_!EV&|&{K`E!U78DLm z8nT?692iDmXBZfy1G&Ghs zz$Q1|)s^R1|7tlHMqiSh1U&Xv5)+9RnEXVvr{lBgn+gqw^WdpE~}^jqI}fmRy>OVIh%K{}?H+w3Xl? z7Ug?BZ^`;!auA0f?j^9}iM7EK52^!UYti;!6Bw2^+(5a38w(H-M8p{BEVhIVc=m$T zNP#<;Ghhhw^S1;^5NrOGF?Atd!+5xvxp~LGP~tYIUhqJhyXI-Mp&4KPzW){+i18#o z83$r@LoxZKuRP}-)Fgcf5Qqw($C(%zH5t`CJ?oQ`r(jxh+G|g4vtAVf6dp)YRay!U zq3P-A=@R$-2m1dH=}&l9wUDh;Q014PBg9g4mJekFfL?go3OtR2Mp)S3(fzxmJgg5M zw6u%I^G`#~5#;AT0}u;MYevQ!>{dfocsz-;k}A|U=q=^tdq7H^@!W2dNXZjpF0JNj z!xM$$)@4{cE`a0C&Ke3j|Nj}4dkeoWN4W&ENQjn#KH0zhxvKQgiyRysnU{0QF)K%O zQ3u$&`nok-{&Lv%B*}XU*lJUiX#cyc12Qs-|2<;j@r(Hhi~B`nD!vDB0V|(=vrs%% z`awrrx@d@v4}}M(91FyUk4%-V>XJ}r^gW&;!Zh!rREa{(+Iw+WJ4_#2@*RV^=$|Y2 z=#8kBY*QEIvlW3f8JHKrk4}IJtO!wcbkxX^qS4J2St?r*vV>$CA{honvP41| zOS&ja(!%$Qet-P@Q!{49JMZ&6pXHqMIfd(hQY5U$1-cYACWpD!;QCxzbVC0!m59%F z>t@|h0m?VFzO4d2mSrse0T=pz>}7HMUb$=kyi_+d@>+9=TCE9~Uo$2S& zj^$my)~mQX)`yQeXzD4)7Ouy4Ns7<2xb>>r=u1(4d!?0$FGlBOi;Gofe%su}H7 znqQG5B+q$w1#U&ncTxbmXO>5Cn7lIUgUS2)Lq#c1Z-3d-hKEr%y)*|i<4^b1wy|&u zlGu=E4UdA#{_pG7;J5o4%S6A-OuF!HIlqThHV|g-;N?O7%|r{X|Mx^Ztt-vn{|c*N zB#dd|6HN14{FJ((Y{1WNmp-0^dWYVGx+PX%akvmw$+o)UG`Ah=)n9hv(Jq^}A{7yA zI62&Q&qIUmPm?VLzU>sO5tOjV+4{Z*1Gqt3C+X{wj{} zw{I;0gvZHSeA5YxnTZa|Mumz3CrEt|K!|hmZ9qh0z{(JLAdX-So_D@_@t!21*6A}= z!>X(aW1LkYX91Xkx7XL_AUFOopTRJ&*{(sQ==7<@ zIbAAZ0%54TIVV28&QO4tk`;6zY|WcX)J$%IYG)arGz zcF7p2s+xV1@`)GFW+>B{Ck`HDAdm3a;A{|5d1u{g#4di9po3hfBc@&D6?36bT9VY8 z+n`uMylF`7_iTTG-hK4t9aj&JcYS@V!3O9-L@Olmd%gLS-JxsA_wL<8Q7WQ*g5zqP zvabSu5o+I7h)M!2_f7kO<^P-^cA=Ewv8(fNKXaxPl?+n*VnN-$hlPb-(PSX1c3jJ} zb#M5o-t!pl6ySVPbVFW^&LCQ4drkPfR@zF=Dk@UZ)lF#_va(+M9Mjo5n(k_eS!5H4 zK&jN#>gmCXA|t0yp5%+vj|s7~wl3;9MX5NivYQWWm#l)q@XP|mlyd#5Y@~krucAXz zF?9MqofM<~M+L&dq7S2@xE138C!9mq2rm&*#ocaU^jw01`U>t%8UD26w|L|s*m>y! zF8?*sWL@3;ii)HS-8#ZNbtj^b($3;^b9nEkZt)U8V(59 zf~(L^DbxM@{9rc%k%du8aP)cxBHfL_IV(-X3 z&;zgyh~01KXCzI0J}?l1N|NMsnfw))!WJp1gmK5~;@1@a#t;?YXMA`k)gzV%!mL{j zLw*4Px2V`3JiNS`rw$fKY~H*XK=tk zYOBUgRj$T-H%gPhO7JE!f@oCTaLXJIf_R30<#_dl9i5>it$W|m(Bu^r#fTH#UzpD` zV`BIt+TaM)SPz39Kz%t^Wu1Wvpn+$euBtBE7s!9C)kEW&zGU6)SrM1kWDJvOd-iR2 z*Edrv8Ty^1<&f5q1YOOZJ?w!tit#DC0u;Ay<^CjzC6$bTApsm_xG>%a5mi=}jEBeP zwl-sIMAfb2F%U5Tx(8YGHWx0dmF8Ve1Db@E^%5r!|EkMC8)`^*#*=OB?I%Y@{6D@} z3A6x)Gs#+k?Pwa!dCg*7Ta)up27s1frrRI44f-wQD>nB#iHE-T|1OTNbAC7K08pAJ zKr1K={s&ZB-!<>r+hg#PlA`O)Y{U*`UC=398>$yGQ&S+*u?@=?1Y0klx;Pw|fP@1B zH8QfWs3<2tUxT4to}NzLDuiJmT?K`_)KmxT5o|30Z%=s5SXC+QTw!UcW_eb3U95@5 zVD<6h+Om^>jsREs5BriQPjX=>htdV3elXl%sy9Eu9e6AirgwIT&J!*U*IT7WBxhq%;+ zO?OR>kN3_qt{7f6#@lnK?olY7k^MXF42^R;I2ypHzilcL~nhv=;{sfrR zhzO22UDe>=U`8eGM(JacQ@FUKq$D31UT(y+z?l)ou=nzLeq>nA=48jj5~F_s0fIG% zXcInu0W+6jhkcFkmuP`+VAxVy8&-b;esdGBpQoU_xz213-hJ9i$btb=IQc6lt}Pl8BHUfyxo;()}Myu!xI50jRIY-L3S=)VW1+0&`nlMEbqkh^UQ z$$ib+c|Lurzf5HnGXzeD*d<=HPXs0Z<@$xd9?FSkjX)x(IpCifLh5qL;w11azld9Hyi0yjZKi4d}xcmk6$z`o7bqsy@_cDbQ1`MO~rX5NNAeneX{S2(Q`b3m#& ze%lnpc)`8s1iA{28!5ou^DN$om4D;oS1^TmjEX_kcOBtmX3D_Cikk1f{{AsbGO;1T z)!Eq$$k`>#t*vuz-r$5D^RK_w2?>om4OfX`oC@|*VaoYgS$afS!R!5_@Q@n?7~4eU z=0pQVHahB7fP2=^FbXfcomQfRom^9)Y~MA(9`5oP;VwCc0oxJ_xwu0a9J{l)V;X|u zY9}W>v(0ndeALy|0nMOs|63RU9XV|*Sl*p0&d=AE7BrVLT6QQ9{WR`tlK>PDa0Kee zFMjZViSGd@hIsZ;)25r~WvdrP?!aBgd;12MSGsq~c~^ogfjYRjc*f742ySzEQ;MEFa%9cT z2JlmLBuxzs_`ji|Gvz7E$k;0*=-E%RQjCYLd0F&9iDJ?Cz}ye>MObbH%b+MDBVK83 zX_;Us-Q3z5;O%|r@ZpJalWp4yu#|{%;t}YFZzm@9baiF=%O`0ZE)V{xK|YD^j`$9` z7!&F`I?1$0u7M<{Fumx+#K@%uV$j!*^++xd457fo;DLFzDpaYmvXR!@L?R#>Rr)zN zn_^CrnVQ*n^ypCtUtUBCc$aJK*s)$hP*~W)#Kh6nbxYhKPz3-*e#8GST3+IjX&@OX zVEcwwNGeKn$`h!QqVG_nqGAt7L723k_!3Y^b+ry+rG(VK2bAb2B8h2GGb+sU5N6r_ z+^R1?hY(yCrI9EZid|kx8ajsw-4KOWy6Qc@yr_bKx z{KaGlOoyR7ny>=O^BF)e?-LL6=)!`=4PX4X9jeN_~pEjYcS1w z`Eot;`?8k}-bkO9nO`Ood83ds)4&3-ZiLFdv!a)<2aDuV(u%&%w!pbvw4 zV*Lq|c~}$1Xl~!|Pm^%Iw#lQTmFr3_$$20e4}U$!oPckTnudmCM3AgeqOpwB8Pwc4 zv8au$?d*VI87eMLH*X^iq3n^2Acl(*B^4A-zpwC1lBDXVm-Yuyt>1yZAe6IC4%(y%nTSr zDxayRhT+J-=#Gzzw&!dnjzMZ0v?fREF&e+Tyc=~N#3R5JX=!Qkk*i`KM#+_If!;@m zdeqKN8_hyQ8(6)kX9gA-TT-IurKP1I=*Dq#r9eM^@803O3G5S;VHk0qpP%RD<%Pfl zSr18pcg-4U`Kic;AxOK8sX#aJs;bz$U#Ec4pE+~p@N0sS`1aKzgrkPiK5&TzoSmTT z#5s2*Mvz2;^{B!0W|k%+IyxHbjuQ>TOwiB~ZQ$h@yn9Pz^(KyW;7de{6Q>Qj@QnsS zmQYy{B_3pJ`S|#hY%!RpxVUIwjvvHZgz=Fh{~VCI{?9+(zVA0THwQZeQz&d2-@;8uOWbCBBG*LJey?8t@dSSKfoUpbFn?~ zF)pv}TA7-L^X-JUsIf60I*rLBCC6|~8@RjYWH$U*SOBfAq@vRP>eVEQLX>o2 z*9qsW(YLg*F-lT981i|b>U@Bs;|C9|HU{GpuzrM5c(-i#8}GpCa;5&1jF{}Ns&soJ lJysq1bZ{eBDuM%HiFd!iki~e`;|4q=9Nc?EKcDOn{yz=4B`^R0 diff --git a/doc/auxiliary-tools.md b/doc/auxiliary-tools.md index 9f6b3d7d..5871b897 100644 --- a/doc/auxiliary-tools.md +++ b/doc/auxiliary-tools.md @@ -80,36 +80,17 @@ With $19\times 19\times 19$ mesh: :width: 25% ``` -###General options - -#### `--pa` - -See {ref}`pa_option`. - -#### `-c` - -Unit cell filename is specified with this option, e.g., `-c POSCAR-unitcell`. - -#### `--qe` - -Let `phono3py-kaccum` read a QE (pw) unit cell file with `-c` option, for -example: - -```bash -% phono3py-kaccum --qe kappa-m191919.hdf5 -``` +That calculated by QE with $19\times 19\times 19$ mesh: ```{image} Si-kaccum-pwscf.png :width: 25% ``` -#### `--crystal` - -Analogous to `--qe`, but to be used with the CRYSTAL interface. +###General options -#### `--turbomole` +#### `--pa` -Analogous to `--qe`, but to be used with the TURBOMOLE interface +See {ref}`pa_option`. #### `--temperature` diff --git a/doc/changelog.md b/doc/changelog.md index e3fb9719..02334d60 100644 --- a/doc/changelog.md +++ b/doc/changelog.md @@ -6,9 +6,9 @@ This is a major version release. There are backward-incompatible changes. -- Method to transform supercell third-order force constants fc3 in real to - reciprocal space was changed as described at Version 2.9.0 changelog below. - This results in the change of results with respect to those obtained by +- Calculation method to transform supercell third-order force constants fc3 in + real to reciprocal space was changed as described at Version 2.9.0 changelog + below. This results in the change of results with respect to those obtained by phono3py version 2. To emulate v2 behaviour, use `--v2` option in phono3py command line script. For `Phono3py` class , `make_r0_average=True` (default) when instantiating it, and similarly for `phono3py.load` function. diff --git a/phono3py/cui/kaccum_script.py b/phono3py/cui/kaccum_script.py index 0937c3e7..f838ac2a 100644 --- a/phono3py/cui/kaccum_script.py +++ b/phono3py/cui/kaccum_script.py @@ -129,14 +129,14 @@ def _assert_grid_in_hdf5( def _get_calculator(args): """Return calculator name.""" interface_mode = None - if args.qe_mode: - interface_mode = "qe" - elif args.crystal_mode: - interface_mode = "crystal" - elif args.abinit_mode: - interface_mode = "abinit" - elif args.turbomole_mode: - interface_mode = "turbomole" + # if args.qe_mode: + # interface_mode = "qe" + # elif args.crystal_mode: + # interface_mode = "crystal" + # elif args.abinit_mode: + # interface_mode = "abinit" + # elif args.turbomole_mode: + # interface_mode = "turbomole" return interface_mode @@ -157,16 +157,6 @@ def _read_files(args): return cell, f -def _read_files_by_collect_cell_info(cell_filename, interface_mode): - cell_info = collect_cell_info( - interface_mode=interface_mode, - cell_filename=cell_filename, - supercell_matrix=np.eye(3, dtype=int), - phonopy_yaml_cls=Phono3pyYaml, - ) - return cell_info - - def _get_mode_property(args, f_kappa): """Read property data from hdf5 file object.""" if args.pqj: @@ -206,14 +196,14 @@ def _get_parser(): default=None, help="Same as PRIMITIVE_AXES tags", ) - parser.add_argument( - "-c", - "--cell", - dest="cell_filename", - metavar="FILE", - default=None, - help="Read unit cell", - ) + # parser.add_argument( + # "-c", + # "--cell", + # dest="cell_filename", + # metavar="FILE", + # default=None, + # help="Read unit cell", + # ) parser.add_argument( "--gv", action="store_true", help="Calculate for gv_x_gv (tensor)" ) @@ -271,24 +261,6 @@ def _get_parser(): action="store_true", help="Use smearing method (only for scalar density)", ) - parser.add_argument( - "--qe", "--pwscf", dest="qe_mode", action="store_true", help="Invoke Pwscf mode" - ) - parser.add_argument( - "--crystal", - dest="crystal_mode", - action="store_true", - help="Invoke CRYSTAL mode", - ) - parser.add_argument( - "--abinit", dest="abinit_mode", action="store_true", help="Invoke Abinit mode" - ) - parser.add_argument( - "--turbomole", - dest="turbomole_mode", - action="store_true", - help="Invoke TURBOMOLE mode", - ) parser.add_argument( "--no-gridsym", dest="no_gridsym", @@ -363,8 +335,12 @@ def main(): 'Use of "phono3py-kaccum CRYSTAL_STRUCTURE_FILE" is not supported.' ) else: - interface_mode = _get_calculator(args) - cell_info = _read_files_by_collect_cell_info(args.cell_filename, interface_mode) + cell_info = collect_cell_info( + supercell_matrix=np.eye(3, dtype=int), + phonopy_yaml_cls=Phono3pyYaml, + ) + cell_filename = cell_info["optional_structure_info"][0] + print(f'# Crystal structure was read from "{cell_filename}".') cell = cell_info["unitcell"] phpy_yaml = cell_info.get("phonopy_yaml", None) if phpy_yaml is not None: diff --git a/phono3py/cui/phono3py_argparse.py b/phono3py/cui/phono3py_argparse.py index 0cbdc115..c15677d2 100644 --- a/phono3py/cui/phono3py_argparse.py +++ b/phono3py/cui/phono3py_argparse.py @@ -282,7 +282,7 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False): help="Read third order force constants", ) parser.add_argument( - "--v3", + "--v2", dest="is_fc3_r0_average", action="store_true", default=False, diff --git a/phono3py/other/kaccum.py b/phono3py/other/kaccum.py index 520cf0a7..edc81f70 100644 --- a/phono3py/other/kaccum.py +++ b/phono3py/other/kaccum.py @@ -209,15 +209,34 @@ def run_prop_dos( mode_prop, ir_grid_map, ir_grid_points, - num_sampling_points, + num_sampling_points: int, bz_grid: BZGrid, ): - """Run DOS-like calculation.""" + """Run DOS-like calculation. + + This is a simple wrapper of KappsDOSTHM. + + Parameters + ---------- + frequencies: + Frequencies at ir-grid points. + mode_prop: + Properties at ir-grid points. + ir_grid_map: + Obtained by get_ir_grid_points(bz_grid)[2]. + ir_grid_points: + Obtained by get_ir_grid_points(bz_grid)[0]. + num_sampling_points: + Number of sampling points in horizontal axis. + bz_grid: + BZ grid. + + """ kappa_dos = KappaDOSTHM( mode_prop, frequencies, bz_grid, - ir_grid_points, + bz_grid.bzg2grg[ir_grid_points], ir_grid_map=ir_grid_map, num_sampling_points=num_sampling_points, ) @@ -227,7 +246,12 @@ def run_prop_dos( def run_mfp_dos( - mean_freepath, mode_prop, ir_grid_map, ir_grid_points, num_sampling_points, bz_grid + mean_freepath, + mode_prop, + ir_grid_map, + ir_grid_points, + num_sampling_points: int, + bz_grid: BZGrid, ): """Run DOS-like calculation for mean free path. @@ -242,7 +266,7 @@ def run_mfp_dos( mode_prop[i : i + 1, :, :], mean_freepath[i], bz_grid, - ir_grid_points, + bz_grid.bzg2grg[ir_grid_points], ir_grid_map=ir_grid_map, num_sampling_points=num_sampling_points, ) diff --git a/test/other/test_kaccum.py b/test/other/test_kaccum.py index a88b40c2..79e53d6c 100644 --- a/test/other/test_kaccum.py +++ b/test/other/test_kaccum.py @@ -7,7 +7,13 @@ import numpy as np from phono3py import Phono3py -from phono3py.other.kaccum import GammaDOSsmearing, KappaDOSTHM, get_mfp +from phono3py.other.kaccum import ( + GammaDOSsmearing, + KappaDOSTHM, + get_mfp, + run_mfp_dos, + run_prop_dos, +) from phono3py.phonon.grid import get_ir_grid_points @@ -291,6 +297,105 @@ def test_GammaDOSsmearing(nacl_pbe: Phono3py): ) +def test_run_prop_dos(si_pbesol: Phono3py): + ph3 = si_pbesol + ph3.mesh_numbers = [7, 7, 7] + ph3.init_phph_interaction() + ph3.run_thermal_conductivity( + temperatures=[ + 300, + ] + ) + bz_grid = ph3.grid + ir_grid_points, _, ir_grid_map = get_ir_grid_points(bz_grid) + tc = ph3.thermal_conductivity + + kdos, sampling_points = run_prop_dos( + tc.frequencies, tc.mode_kappa[0], ir_grid_map, ir_grid_points, 10, bz_grid + ) + mean_freepath = get_mfp(tc.gamma[0], tc.group_velocities) + mfp, sampling_points_mfp = run_mfp_dos( + mean_freepath, tc.mode_kappa[0], ir_grid_map, ir_grid_points, 10, bz_grid + ) + + # print(",".join([f"{v:10.5f}" for v in kdos[0, :, :, 0].ravel()])) + ref_kdos = [ + 0.00000, + 0.00000, + 2.19162, + 5.16697, + 28.22125, + 18.97280, + 58.56343, + 12.19206, + 69.05896, + 3.47035, + 73.17626, + 1.48915, + 74.74544, + 0.43485, + 75.87064, + 1.74135, + 79.08179, + 2.30428, + 81.21678, + 0.00000, + ] + # print(",".join([f"{v:10.5f}" for v in mfp[0, :, :, 0].ravel()])) + ref_mfp = [ + 0.00000, + 0.00000, + 29.19150, + 0.02604, + 42.80717, + 0.01202, + 52.09457, + 0.01158, + 61.79908, + 0.01140, + 69.49177, + 0.00784, + 74.57499, + 0.00501, + 77.99145, + 0.00364, + 80.33477, + 0.00210, + 81.21678, + 0.00000, + ] + # print(",".join([f"{v:10.5f}" for v in sampling_points[0]])) + ref_sampling_points = [ + -0.00000, + 1.69664, + 3.39328, + 5.08992, + 6.78656, + 8.48320, + 10.17984, + 11.87648, + 13.57312, + 15.26976, + ] + # print(",".join([f"{v:10.5f}" for v in sampling_points_mfp[0]])) + ref_sampling_points_mfp = [ + 0.00000, + 803.91710, + 1607.83420, + 2411.75130, + 3215.66841, + 4019.58551, + 4823.50261, + 5627.41971, + 6431.33681, + 7235.25391, + ] + np.testing.assert_allclose(ref_kdos, kdos[0, :, :, 0].ravel(), atol=1e-2) + np.testing.assert_allclose(ref_mfp, mfp[0, :, :, 0].ravel(), atol=1e-2) + np.testing.assert_allclose(ref_sampling_points, sampling_points[0], atol=1e-4) + np.testing.assert_allclose(ref_sampling_points_mfp, sampling_points_mfp[0], rtol=10) + + def _calculate_kappados( ph3: Phono3py, mode_prop: np.ndarray, From e50cecc1cf8bcd4957fcff29539d04ae65b56225 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Fri, 19 Apr 2024 16:51:55 +0900 Subject: [PATCH 15/16] Loosen testing criteria for test_run_prop_dos --- doc/auxiliary-tools.md | 12 +++---- phono3py/cui/kaccum_script.py | 67 ----------------------------------- test/other/test_kaccum.py | 4 +-- 3 files changed, 7 insertions(+), 76 deletions(-) diff --git a/doc/auxiliary-tools.md b/doc/auxiliary-tools.md index 5871b897..c54c9741 100644 --- a/doc/auxiliary-tools.md +++ b/doc/auxiliary-tools.md @@ -35,7 +35,7 @@ $$ ### How to use `phono3py-kaccum` -Let's computer lattice thermal conductivity of Si using the `Si-PBEsol` example +Let's compute lattice thermal conductivity of Si using the `Si-PBEsol` example found in the example directory. ```bash @@ -50,7 +50,9 @@ follows: ``` `kappa-m111111.hdf5` is the output file of thermal conductivity calculation, -which is passed to `phono3py-kaccum` as the first argument. +which is passed to `phono3py-kaccum` as the first argument. `phono3py_disp.yaml` +or `phono3py.yaml` is necessary to be located under the current directory to +read the crystal structure. The format of the output is as follows: The first column gives frequency in THz, and the second to seventh columns give the cumulative lattice thermal @@ -86,11 +88,7 @@ That calculated by QE with $19\times 19\times 19$ mesh: :width: 25% ``` -###General options - -#### `--pa` - -See {ref}`pa_option`. +### General options #### `--temperature` diff --git a/phono3py/cui/kaccum_script.py b/phono3py/cui/kaccum_script.py index f838ac2a..149bf244 100644 --- a/phono3py/cui/kaccum_script.py +++ b/phono3py/cui/kaccum_script.py @@ -6,13 +6,7 @@ import h5py import numpy as np from phonopy.cui.collect_cell_info import collect_cell_info -from phonopy.cui.settings import fracval from phonopy.interface.calculator import read_crystal_structure -from phonopy.structure.cells import ( - get_primitive, - get_primitive_matrix, - guess_primitive_matrix, -) from phonopy.structure.symmetry import Symmetry from phono3py.interface.phono3py_yaml import Phono3pyYaml @@ -187,23 +181,6 @@ def _get_mode_property(args, f_kappa): def _get_parser(): """Return args of ArgumentParser.""" parser = argparse.ArgumentParser(description="Show unit cell volume") - parser.add_argument( - "--pa", - "--primitive-axis", - "--primitive-axes", - nargs="+", - dest="primitive_matrix", - default=None, - help="Same as PRIMITIVE_AXES tags", - ) - # parser.add_argument( - # "-c", - # "--cell", - # dest="cell_filename", - # metavar="FILE", - # default=None, - # help="Read unit cell", - # ) parser.add_argument( "--gv", action="store_true", help="Calculate for gv_x_gv (tensor)" ) @@ -272,34 +249,6 @@ def _get_parser(): return args -def _analyze_primitive_matrix_option(args, unitcell=None): - """Analyze --pa option argument.""" - if args.primitive_matrix is not None: - if type(args.primitive_matrix) is list: - _primitive_matrix = " ".join(args.primitive_matrix) - else: - _primitive_matrix = args.primitive_matrix.strip() - - if _primitive_matrix.lower() == "auto": - primitive_matrix = "auto" - elif _primitive_matrix.upper() in ("P", "F", "I", "A", "C", "R"): - primitive_matrix = _primitive_matrix.upper() - elif len(_primitive_matrix.split()) != 9: - raise SyntaxError("Number of elements in --pa option argument has to be 9.") - else: - primitive_matrix = np.reshape( - [fracval(x) for x in _primitive_matrix.split()], (3, 3) - ) - if np.linalg.det(primitive_matrix) < 1e-8: - raise SyntaxError("Primitive matrix has to have positive determinant.") - - pmat = get_primitive_matrix(primitive_matrix) - if unitcell is not None and isinstance(pmat, str) and pmat == "auto": - return guess_primitive_matrix(unitcell) - else: - return pmat - - def main(): """Calculate kappa spectrum. @@ -310,12 +259,6 @@ def main(): % phono3py-kaccum kappa-m111111.hdf5 ``` - Old style usage - --------------- - ``` - % phono3py-kaccum --pa="F" -c POSCAR-unitcell kappa-m111111.hdf5 |tee kaccum.dat - ``` - Plot by gnuplot --------------- ``` @@ -324,9 +267,6 @@ def main(): gnuplot> p "kaccum.dat" i 30 u 1:2 w l, "kaccum.dat" i 30 u 1:8 w l ``` - With phono3py.yaml type file as crystal structure, primitive matrix is - unnecessary to set. - """ args = _get_parser() primitive = None @@ -347,14 +287,7 @@ def main(): primitive = cell_info["phonopy_yaml"].primitive if primitive is None: primitive = cell - else: - primitive_matrix = _analyze_primitive_matrix_option(args, unitcell=cell) f_kappa = h5py.File(args.filenames[0], "r") - if primitive is None: - if primitive_matrix is None: - primitive = cell - else: - primitive = get_primitive(cell, primitive_matrix) if "grid_matrix" in f_kappa: mesh = np.array(f_kappa["grid_matrix"][:], dtype="int_") diff --git a/test/other/test_kaccum.py b/test/other/test_kaccum.py index 79e53d6c..c6bd1eac 100644 --- a/test/other/test_kaccum.py +++ b/test/other/test_kaccum.py @@ -390,8 +390,8 @@ def test_run_prop_dos(si_pbesol: Phono3py): 6431.33681, 7235.25391, ] - np.testing.assert_allclose(ref_kdos, kdos[0, :, :, 0].ravel(), atol=1e-2) - np.testing.assert_allclose(ref_mfp, mfp[0, :, :, 0].ravel(), atol=1e-2) + np.testing.assert_allclose(ref_kdos, kdos[0, :, :, 0].ravel(), atol=1e-1) + np.testing.assert_allclose(ref_mfp, mfp[0, :, :, 0].ravel(), atol=1e-1) np.testing.assert_allclose(ref_sampling_points, sampling_points[0], atol=1e-4) np.testing.assert_allclose(ref_sampling_points_mfp, sampling_points_mfp[0], rtol=10) From cc7fbfefa59b97255703c2c75b8fb5b753cf76e9 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Fri, 19 Apr 2024 16:59:52 +0900 Subject: [PATCH 16/16] Loosen testing criteria for test_run_prop_dos --- test/other/test_kaccum.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/other/test_kaccum.py b/test/other/test_kaccum.py index c6bd1eac..503d78af 100644 --- a/test/other/test_kaccum.py +++ b/test/other/test_kaccum.py @@ -390,8 +390,8 @@ def test_run_prop_dos(si_pbesol: Phono3py): 6431.33681, 7235.25391, ] - np.testing.assert_allclose(ref_kdos, kdos[0, :, :, 0].ravel(), atol=1e-1) - np.testing.assert_allclose(ref_mfp, mfp[0, :, :, 0].ravel(), atol=1e-1) + np.testing.assert_allclose(ref_kdos, kdos[0, :, :, 0].ravel(), atol=0.5) + np.testing.assert_allclose(ref_mfp, mfp[0, :, :, 0].ravel(), atol=0.5) np.testing.assert_allclose(ref_sampling_points, sampling_points[0], atol=1e-4) np.testing.assert_allclose(ref_sampling_points_mfp, sampling_points_mfp[0], rtol=10)