diff --git a/phono3py/file_IO.py b/phono3py/file_IO.py index 2be4be36..4bd764fc 100644 --- a/phono3py/file_IO.py +++ b/phono3py/file_IO.py @@ -453,7 +453,12 @@ def write_grid_address(grid_address, mesh, filename=None): def write_grid_address_to_hdf5( - grid_address, mesh, grid_mapping_table, compression="gzip", filename=None + grid_address, + mesh, + grid_mapping_table, + bz_grid=None, + compression="gzip", + filename=None, ): """Write grid addresses to grid_address.hdf5.""" suffix = _get_filename_suffix(mesh, filename=filename) @@ -461,6 +466,10 @@ def write_grid_address_to_hdf5( with h5py.File(full_filename, "w") as w: w.create_dataset("version", data=np.string_(__version__)) w.create_dataset("mesh", data=mesh) + if bz_grid is not None and bz_grid.grid_matrix is not None: + w.create_dataset("grid_matrix", data=bz_grid.grid_matrix) + w.create_dataset("P_matrix", data=bz_grid.P) + w.create_dataset("Q_matrix", data=bz_grid.Q) w.create_dataset("grid_address", data=grid_address, compression=compression) w.create_dataset( "grid_mapping_table", data=grid_mapping_table, compression=compression @@ -618,6 +627,7 @@ def write_real_self_energy_to_hdf5( deltas, mesh, epsilon, + bz_grid=None, frequency_points=None, frequencies=None, filename=None, @@ -649,6 +659,10 @@ def write_real_self_energy_to_hdf5( w.create_dataset("version", data=np.string_(__version__)) w.create_dataset("grid_point", data=grid_point) w.create_dataset("mesh", data=mesh) + if bz_grid is not None and bz_grid.grid_matrix is not None: + w.create_dataset("grid_matrix", data=bz_grid.grid_matrix) + w.create_dataset("P_matrix", data=bz_grid.P) + w.create_dataset("Q_matrix", data=bz_grid.Q) w.create_dataset("band_index", data=_band_indices) w.create_dataset("delta", data=deltas) w.create_dataset("temperature", data=temperatures) @@ -700,6 +714,7 @@ def write_spectral_function_to_hdf5( shifts, half_linewidths, mesh, + bz_grid=None, sigma=None, frequency_points=None, frequencies=None, @@ -724,6 +739,10 @@ def write_spectral_function_to_hdf5( w.create_dataset("version", data=np.string_(__version__)) w.create_dataset("grid_point", data=grid_point) w.create_dataset("mesh", data=mesh) + if bz_grid is not None and bz_grid.grid_matrix is not None: + w.create_dataset("grid_matrix", data=bz_grid.grid_matrix) + w.create_dataset("P_matrix", data=bz_grid.P) + w.create_dataset("Q_matrix", data=bz_grid.Q) w.create_dataset("band_index", data=_band_indices) w.create_dataset("spectral_function", data=spectral_functions) w.create_dataset("shift", data=shifts) @@ -884,6 +903,7 @@ def write_collision_eigenvalues_to_hdf5( def write_kappa_to_hdf5( temperature, mesh, + bz_grid=None, frequency=None, group_velocity=None, gv_by_gv=None, @@ -928,7 +948,10 @@ def write_kappa_to_hdf5( w.create_dataset("version", data=np.string_(__version__)) w.create_dataset("temperature", data=temperature) w.create_dataset("mesh", data=mesh) - + if bz_grid is not None and bz_grid.grid_matrix is not None: + w.create_dataset("grid_matrix", data=bz_grid.grid_matrix) + w.create_dataset("P_matrix", data=bz_grid.P) + w.create_dataset("Q_matrix", data=bz_grid.Q) if frequency is not None: if isinstance(frequency, np.floating): w.create_dataset("frequency", data=frequency) @@ -1299,6 +1322,7 @@ def read_pp_from_hdf5( def write_gamma_detail_to_hdf5( temperature, mesh, + bz_grid=None, gamma_detail=None, grid_point=None, triplet=None, @@ -1332,6 +1356,10 @@ def write_gamma_detail_to_hdf5( w.create_dataset("version", data=np.string_(__version__)) w.create_dataset("temperature", data=temperature) w.create_dataset("mesh", data=mesh) + if bz_grid is not None and bz_grid.grid_matrix is not None: + w.create_dataset("grid_matrix", data=bz_grid.grid_matrix) + w.create_dataset("P_matrix", data=bz_grid.P) + w.create_dataset("Q_matrix", data=bz_grid.Q) if gamma_detail is not None: w.create_dataset("gamma_detail", data=gamma_detail, compression=compression) if triplet is not None: @@ -1384,7 +1412,13 @@ def write_gamma_detail_to_hdf5( def write_phonon_to_hdf5( - frequency, eigenvector, grid_address, mesh, compression="gzip", filename=None + frequency, + eigenvector, + grid_address, + mesh, + bz_grid=None, + compression="gzip", + filename=None, ): """Write phonon on grid in its hdf5 file.""" suffix = _get_filename_suffix(mesh, filename=filename) @@ -1393,6 +1427,10 @@ def write_phonon_to_hdf5( with h5py.File(full_filename, "w") as w: w.create_dataset("version", data=np.string_(__version__)) w.create_dataset("mesh", data=mesh) + if bz_grid is not None and bz_grid.grid_matrix is not None: + w.create_dataset("grid_matrix", data=bz_grid.grid_matrix) + w.create_dataset("P_matrix", data=bz_grid.P) + w.create_dataset("Q_matrix", data=bz_grid.Q) w.create_dataset("grid_address", data=grid_address, compression=compression) w.create_dataset("frequency", data=frequency, compression=compression) w.create_dataset("eigenvector", data=eigenvector, compression=compression) diff --git a/phono3py/phonon3/conductivity.py b/phono3py/phonon3/conductivity.py index e885b0db..17e5f5a1 100644 --- a/phono3py/phonon3/conductivity.py +++ b/phono3py/phonon3/conductivity.py @@ -278,6 +278,11 @@ class Conductivity(ConductivityBase): ) return self.mesh_numbers + @property + def bz_grid(self): + """Return GR-grid.""" + return self._pp.bz_grid + @property def mode_heat_capacities(self): """Return mode heat capacity at constant volume at grid points. diff --git a/phono3py/phonon3/conductivity_LBTE.py b/phono3py/phonon3/conductivity_LBTE.py index 33c93454..87f5576f 100644 --- a/phono3py/phonon3/conductivity_LBTE.py +++ b/phono3py/phonon3/conductivity_LBTE.py @@ -1486,6 +1486,7 @@ def _write_kappa( sigmas = lbte.sigmas sigma_cutoff = lbte.sigma_cutoff_width mesh = lbte.mesh_numbers + bz_grid = lbte.bz_grid weights = lbte.grid_weights frequencies = lbte.frequencies ave_pp = lbte.averaged_pp_interaction @@ -1519,6 +1520,7 @@ def _write_kappa( write_kappa_to_hdf5( temperatures, mesh, + bz_grid=bz_grid, frequency=frequencies, group_velocity=gv, gv_by_gv=gv_by_gv, diff --git a/phono3py/phonon3/conductivity_RTA.py b/phono3py/phonon3/conductivity_RTA.py index ae946410..989b5f50 100644 --- a/phono3py/phonon3/conductivity_RTA.py +++ b/phono3py/phonon3/conductivity_RTA.py @@ -35,6 +35,7 @@ import sys +from typing import Optional import numpy as np from phonopy.phonon.group_velocity import GroupVelocity from phonopy.structure.tetrahedron_method import TetrahedronMethod @@ -50,7 +51,7 @@ from phono3py.phonon3.conductivity import write_pp as _write_pp from phono3py.phonon3.imag_self_energy import ImagSelfEnergy, average_by_degeneracy from phono3py.phonon3.interaction import Interaction from phono3py.phonon3.triplets import get_all_triplets -from phono3py.phonon.grid import get_grid_points_by_rotations +from phono3py.phonon.grid import _get_ir_grid_map, get_grid_points_by_rotations class ConductivityRTA(Conductivity): @@ -676,15 +677,21 @@ def get_thermal_conductivity_RTA( def _write_gamma_detail( - br, interaction, i, compression="gzip", filename=None, verbose=True + br: ConductivityRTA, + interaction: Interaction, + i: int, + compression: str = "gzip", + filename: Optional[str] = None, + verbose: bool = True, ): gamma_detail = br.get_gamma_detail_at_q() - temperatures = br.get_temperatures() - mesh = br.get_mesh_numbers() - grid_points = br.get_grid_points() + temperatures = br.temperatures + mesh = br.mesh_numbers + bz_grid = br.bz_grid + grid_points = br.grid_points gp = grid_points[i] - sigmas = br.get_sigmas() - sigma_cutoff = br.get_sigma_cutoff_width() + sigmas = br.sigmas + sigma_cutoff = br.sigma_cutoff_width triplets, weights, _, _ = interaction.get_triplets_at_q() all_triplets = get_all_triplets(gp, interaction.bz_grid) @@ -693,6 +700,7 @@ def _write_gamma_detail( write_gamma_detail_to_hdf5( temperatures, mesh, + bz_grid=bz_grid, gamma_detail=gamma_detail, grid_point=gp, triplet=triplets, @@ -710,6 +718,7 @@ def _write_gamma_detail( write_gamma_detail_to_hdf5( temperatures, mesh, + bz_grid=bz_grid, gamma_detail=gamma_detail[:, :, k, :, :], grid_point=gp, triplet=triplets, @@ -723,7 +732,14 @@ def _write_gamma_detail( ) -def _write_gamma(br, interaction, i, compression="gzip", filename=None, verbose=True): +def _write_gamma( + br: ConductivityRTA, + interaction: Interaction, + i: int, + compression: str = "gzip", + filename: Optional[str] = None, + verbose: bool = True, +): """Write mode kappa related properties into a hdf5 file.""" grid_points = br.grid_points group_velocities = br.group_velocities @@ -731,6 +747,7 @@ def _write_gamma(br, interaction, i, compression="gzip", filename=None, verbose= mode_heat_capacities = br.mode_heat_capacities ave_pp = br.averaged_pp_interaction mesh = br.mesh_numbers + bz_grid = br.bz_grid temperatures = br.temperatures gamma = br.gamma gamma_isotope = br.gamma_isotope @@ -763,6 +780,7 @@ def _write_gamma(br, interaction, i, compression="gzip", filename=None, verbose= write_kappa_to_hdf5( temperatures, mesh, + bz_grid=bz_grid, frequency=frequencies, group_velocity=group_velocities[i], gv_by_gv=gv_by_gv[i], @@ -803,6 +821,7 @@ def _write_gamma(br, interaction, i, compression="gzip", filename=None, verbose= write_kappa_to_hdf5( temperatures, mesh, + bz_grid=bz_grid, frequency=frequencies, group_velocity=group_velocities[i, k], gv_by_gv=gv_by_gv[i, k], @@ -861,7 +880,13 @@ def _show_kappa(br, log_level): print("") -def _write_kappa(br, volume, compression="gzip", filename=None, log_level=0): +def _write_kappa( + br: ConductivityRTA, + volume: float, + compression: str = "gzip", + filename: Optional[str] = None, + log_level: int = 0, +): temperatures = br.temperatures sigmas = br.sigmas sigma_cutoff = br.sigma_cutoff_width @@ -869,6 +894,7 @@ def _write_kappa(br, volume, compression="gzip", filename=None, log_level=0): gamma_isotope = br.gamma_isotope gamma_N, gamma_U = br.get_gamma_N_U() mesh = br.mesh_numbers + bz_grid = br.bz_grid frequencies = br.frequencies gv = br.group_velocities gv_by_gv = br.gv_by_gv @@ -897,6 +923,7 @@ def _write_kappa(br, volume, compression="gzip", filename=None, log_level=0): write_kappa_to_hdf5( temperatures, mesh, + bz_grid=bz_grid, frequency=frequencies, group_velocity=gv, gv_by_gv=gv_by_gv, diff --git a/phono3py/phonon3/imag_self_energy.py b/phono3py/phonon3/imag_self_energy.py index ade96881..92157714 100644 --- a/phono3py/phonon3/imag_self_energy.py +++ b/phono3py/phonon3/imag_self_energy.py @@ -260,7 +260,7 @@ def _get_imag_self_energy_at_gp( _num_frequency_points, scattering_event_class, num_points_in_batch, - interaction, + interaction: Interaction, ise, write_gamma_detail, return_gamma_detail, @@ -270,6 +270,7 @@ def _get_imag_self_energy_at_gp( num_band0 = len(interaction.band_indices) frequencies = interaction.get_phonons()[0] mesh = interaction.mesh_numbers + bz_grid = interaction.bz_grid if write_gamma_detail or return_gamma_detail: triplets, weights, _, _ = interaction.get_triplets_at_q() @@ -329,6 +330,7 @@ def _get_imag_self_energy_at_gp( full_filename = write_gamma_detail_to_hdf5( temperatures, mesh, + bz_grid=bz_grid, gamma_detail=detailed_gamma_at_gp[j], grid_point=gp, triplet=triplets, diff --git a/phono3py/phonon3/real_self_energy.py b/phono3py/phonon3/real_self_energy.py index 4676d9ac..d069a984 100644 --- a/phono3py/phonon3/real_self_energy.py +++ b/phono3py/phonon3/real_self_energy.py @@ -45,10 +45,11 @@ from phono3py.file_IO import ( ) from phono3py.phonon3.imag_self_energy import get_frequency_points from phono3py.phonon.func import bose_einstein +from phono3py.phonon3.interaction import Interaction def get_real_self_energy( - interaction, + interaction: Interaction, grid_points, temperatures, epsilons=None, @@ -134,6 +135,7 @@ def get_real_self_energy( fst = RealSelfEnergy(interaction) mesh = interaction.mesh_numbers + bz_grid = interaction.bz_grid frequencies = interaction.get_phonons()[0] max_phonon_freq = np.amax(frequencies) band_indices = interaction.band_indices @@ -228,6 +230,7 @@ def get_real_self_energy( all_deltas[i, :, j], mesh, fst.epsilon, + bz_grid=bz_grid, frequency_points=_frequency_points, frequencies=frequencies, filename=output_filename, diff --git a/phono3py/phonon3/spectral_function.py b/phono3py/phonon3/spectral_function.py index 83d31953..b7b4ffa9 100644 --- a/phono3py/phonon3/spectral_function.py +++ b/phono3py/phonon3/spectral_function.py @@ -46,11 +46,12 @@ from phono3py.phonon3.imag_self_energy import ( get_frequency_points, run_ise_at_frequency_points_batch, ) +from phono3py.phonon3.interaction import Interaction from phono3py.phonon3.real_self_energy import imag_to_real def run_spectral_function( - interaction, + interaction: Interaction, grid_points, temperatures=None, sigmas=None, @@ -157,6 +158,7 @@ def run_spectral_function( spf.shifts[sigma_i, :, i], spf.half_linewidths[sigma_i, :, i], interaction.mesh_numbers, + interaction.bz_grid, sigma=sigma, frequency_points=spf.frequency_points, frequencies=frequencies[gp], diff --git a/test/conftest.py b/test/conftest.py index ab11dd8f..baa7eb68 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -48,6 +48,28 @@ def si_pbesol(request): ) +@pytest.fixture(scope="session") +def si_pbesol_grg(request): + """Return Phono3py instance of Si 2x2x2. + + * with symmetry + * full fc + * GR-grid + + """ + yaml_filename = os.path.join(current_dir, "phono3py_si_pbesol.yaml") + forces_fc3_filename = os.path.join(current_dir, "FORCES_FC3_si_pbesol") + enable_v2 = request.config.getoption("--v1") + return phono3py.load( + yaml_filename, + forces_fc3_filename=forces_fc3_filename, + store_dense_gp_map=enable_v2, + store_dense_svecs=enable_v2, + use_grg=True, + log_level=1, + ) + + @pytest.fixture(scope="session") def si_pbesol_nosym(request): """Return Phono3py instance of Si 2x2x2. diff --git a/test/phonon3/test_kappa_RTA.py b/test/phonon3/test_kappa_RTA.py index 11de97cf..e170515b 100644 --- a/test/phonon3/test_kappa_RTA.py +++ b/test/phonon3/test_kappa_RTA.py @@ -1,5 +1,6 @@ """Test for Conductivity_RTA.py.""" import numpy as np +from phono3py import Phono3py si_pbesol_kappa_RTA = [107.991, 107.991, 107.991, 0, 0, 0] si_pbesol_kappa_RTA_with_sigmas = [109.6985, 109.6985, 109.6985, 0, 0, 0] @@ -14,6 +15,7 @@ si_pbesol_kappa_RTA_si_nosym = [ 0.283, ] si_pbesol_kappa_RTA_si_nomeshsym = [38.90918, 38.90918, 38.90918, 0, 0, 0] +si_pbesol_kappa_RTA_grg = [37.3176477, 37.3176477, 37.3176477, 0, 0, 0] nacl_pbe_kappa_RTA = [7.72798252, 7.72798252, 7.72798252, 0, 0, 0] nacl_pbe_kappa_RTA_with_sigma = [7.71913708, 7.71913708, 7.71913708, 0, 0, 0] @@ -95,6 +97,35 @@ def test_kappa_RTA_si_nomeshsym(si_pbesol, si_pbesol_nomeshsym): np.testing.assert_allclose(kappa_ref, kappa, atol=0.5) +def test_kappa_RTA_si_grg(si_pbesol_grg: Phono3py): + """Test RTA by Si with GR-grid.""" + mesh = 10 + ph3 = si_pbesol_grg + ph3.mesh_numbers = mesh + ph3.init_phph_interaction() + ph3.run_thermal_conductivity( + temperatures=[ + 300, + ], + ) + kappa = ph3.thermal_conductivity.kappa.ravel() + np.testing.assert_equal( + ph3.thermal_conductivity.bz_grid.grid_matrix, + [[-2, 2, 2], [2, -2, 2], [2, 2, -2]], + ) + np.testing.assert_equal( + ph3.grid.grid_matrix, + [[-2, 2, 2], [2, -2, 2], [2, 2, -2]], + ) + A = ph3.grid.grid_matrix + D_diag = ph3.grid.D_diag + P = ph3.grid.P + Q = ph3.grid.Q + np.testing.assert_equal(np.dot(P, np.dot(A, Q)), np.diag(D_diag)) + + np.testing.assert_allclose(si_pbesol_kappa_RTA_grg, kappa, atol=0.5) + + def test_kappa_RTA_si_N_U(si_pbesol): """Test RTA with N and U scatterings by Si.""" ph3 = si_pbesol