From 6cf1ff84f2bc1a61c0bff9d3dadc23f5530630a1 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Tue, 7 Apr 2026 11:10:55 +0900 Subject: [PATCH 1/4] Refactoring of conductivity module --- phono3py/conductivity/grid_point_data.py | 2 +- .../conductivity/heat_capacity_providers.py | 4 + phono3py/conductivity/lbte_calculator.py | 70 +++++----------- .../ms_smm19/velocity_providers.py | 7 +- phono3py/conductivity/protocols.py | 38 ++++++--- phono3py/conductivity/rta_calculator.py | 81 +++++++------------ phono3py/conductivity/rta_output.py | 2 +- phono3py/conductivity/velocity_providers.py | 6 ++ .../test_phono3py_load_script_ms_smm19.py | 21 ----- 9 files changed, 89 insertions(+), 142 deletions(-) diff --git a/phono3py/conductivity/grid_point_data.py b/phono3py/conductivity/grid_point_data.py index d9cafbba..62de8f5c 100644 --- a/phono3py/conductivity/grid_point_data.py +++ b/phono3py/conductivity/grid_point_data.py @@ -72,7 +72,7 @@ class HeatCapacityResult: """ - heat_capacities: NDArray[np.double] | None = None + heat_capacities: NDArray[np.double] heat_capacity_matrix: NDArray[np.double] | None = None extra: dict[str, Any] = field(default_factory=dict) diff --git a/phono3py/conductivity/heat_capacity_providers.py b/phono3py/conductivity/heat_capacity_providers.py index 4be3c426..e6011bc3 100644 --- a/phono3py/conductivity/heat_capacity_providers.py +++ b/phono3py/conductivity/heat_capacity_providers.py @@ -61,6 +61,8 @@ class ModeHeatCapacityProvider: (``phonon_all_done == True``) before calling ``compute``. """ + produces_heat_capacity_matrix: bool = False + def __init__(self, pp: Interaction): """Init method.""" self._pp = pp @@ -110,6 +112,8 @@ class HeatCapacityMatrixProvider: """ + produces_heat_capacity_matrix: bool = True + def __init__(self, pp: Interaction): """Init method.""" self._pp = pp diff --git a/phono3py/conductivity/lbte_calculator.py b/phono3py/conductivity/lbte_calculator.py index fac0e56e..f4e138a7 100644 --- a/phono3py/conductivity/lbte_calculator.py +++ b/phono3py/conductivity/lbte_calculator.py @@ -145,7 +145,6 @@ def __init__( self._gamma_boundary: NDArray[np.double] | None = None self._vm_by_vm: NDArray[np.cdouble] | None = None self._heat_capacity_matrix: NDArray[np.double] | None = None - self._extra: dict[str, Any] = {} # Allocate arrays. self._allocate_values() @@ -499,12 +498,29 @@ def _allocate_values(self) -> None: num_temp = len(self._context.temperatures) num_gp = len(self._context.ir_grid_points) num_band0 = len(self._context.band_indices) + num_band = self._context.frequencies.shape[1] self._gamma = np.zeros( (num_sigma, num_temp, num_gp, num_band0), order="C", dtype="double" ) self._gv = np.zeros((num_gp, num_band0, 3), order="C", dtype="double") self._cv = np.zeros((num_temp, num_gp, num_band0), order="C", dtype="double") + if self._isotope_provider is not None: + self._gamma_iso = np.zeros( + (num_sigma, num_gp, num_band0), order="C", dtype="double" + ) + if self._is_full_pp: + self._averaged_pp_interaction = np.zeros( + (num_gp, num_band0), order="C", dtype="double" + ) + if self._velocity_provider.produces_vm_by_vm: + self._vm_by_vm = np.zeros( + (num_gp, num_band0, num_band, 6), order="C", dtype="complex128" + ) + if self._cv_provider.produces_heat_capacity_matrix: + self._heat_capacity_matrix = np.zeros( + (num_temp, num_gp, num_band0, num_band), order="C", dtype="double" + ) def _build_grid_point_aggregates(self) -> GridPointAggregates: """Build GridPointAggregates for accumulator.finalize().""" @@ -518,7 +534,6 @@ def _build_grid_point_aggregates(self) -> GridPointAggregates: gamma_elph=self._gamma_elph, vm_by_vm=self._vm_by_vm, heat_capacity_matrix=self._heat_capacity_matrix, - extra=self._extra, ) def _run_at_grid_point(self, i_gp: int) -> None: @@ -532,72 +547,29 @@ def _run_at_grid_point(self, i_gp: int) -> None: # Mode heat capacities. cv_result = self._cv_provider.compute(gp_input, self._context.temperatures) - cv = cv_result.heat_capacities # (num_temp, num_band0) # Collision matrix row + gamma (Stage 1). collision_result = self._collision_provider.compute(gp_input) # Store per-grid-point data in calculator. - assert self._gamma is not None - assert self._gv is not None - assert self._cv is not None self._gamma[:, :, i_gp, :] = collision_result.gamma self._gv[i_gp] = gv - self._cv[:, i_gp, :] = cv + self._cv[:, i_gp, :] = cv_result.heat_capacities # Isotope scattering (optional). if self._isotope_provider is not None: gamma_iso = self._isotope_provider.compute(gp_input) - gamma_iso = gamma_iso[:, self._context.band_indices] - if self._gamma_iso is None: - self._gamma_iso = np.zeros( - (len(self._context.sigmas),) - + (len(self._context.ir_grid_points),) - + gamma_iso.shape[1:], - order="C", - dtype="double", - ) - self._gamma_iso[:, i_gp, :] = gamma_iso + self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._context.band_indices] # Averaged pp interaction (optional). if self._is_full_pp and collision_result.averaged_pp is not None: - ave_pp = collision_result.averaged_pp - if self._averaged_pp_interaction is None: - self._averaged_pp_interaction = np.zeros( - (len(self._context.ir_grid_points),) + ave_pp.shape, - order="C", - dtype="double", - ) - self._averaged_pp_interaction[i_gp] = ave_pp + self._averaged_pp_interaction[i_gp] = collision_result.averaged_pp # Store vm_by_vm and heat_capacity_matrix (for Wigner/Kubo). if vel_result.vm_by_vm is not None: - if self._vm_by_vm is None: - num_ir = len(self._context.ir_grid_points) - self._vm_by_vm = np.zeros( - (num_ir,) + vel_result.vm_by_vm.shape, dtype="complex128", order="C" - ) self._vm_by_vm[i_gp] = vel_result.vm_by_vm if cv_result.heat_capacity_matrix is not None: - if self._heat_capacity_matrix is None: - shape = cv_result.heat_capacity_matrix.shape - num_ir = len(self._context.ir_grid_points) - self._heat_capacity_matrix = np.zeros( - (shape[0], num_ir) + shape[1:], - dtype="double", - order="C", - ) - self._heat_capacity_matrix[:, i_gp, :, :] = cv_result.heat_capacity_matrix - - # Store velocity extra data (e.g. velocity_operator for Wigner). - if vel_result.extra: - for key, val in vel_result.extra.items(): - if key not in self._extra: - num_ir = len(self._context.ir_grid_points) - self._extra[key] = np.zeros( - (num_ir,) + val.shape, dtype=val.dtype, order="C" - ) - self._extra[key][i_gp] = val + self._heat_capacity_matrix[:, i_gp] = cv_result.heat_capacity_matrix # Accumulate collision matrix row. self._accumulator.store(i_gp, collision_result) diff --git a/phono3py/conductivity/ms_smm19/velocity_providers.py b/phono3py/conductivity/ms_smm19/velocity_providers.py index 77fd2205..f0e211fc 100644 --- a/phono3py/conductivity/ms_smm19/velocity_providers.py +++ b/phono3py/conductivity/ms_smm19/velocity_providers.py @@ -56,6 +56,9 @@ class VelocityOperatorProvider: """ + produces_gv_by_gv: bool = False + produces_vm_by_vm: bool = True + def __init__( self, pp: Interaction, @@ -109,14 +112,10 @@ def compute(self, gp: GridPointInput) -> VelocityResult: gv = self._get_group_velocities(gp, gv_op_full) vm_by_vm, kstar_order = self._get_gv_by_gv_operator(gp, gv_op) - num_band0 = vm_by_vm.shape[0] - gv_by_gv = np.real(vm_by_vm[np.arange(num_band0), np.arange(num_band0)]) return VelocityResult( group_velocities=gv, - gv_by_gv=gv_by_gv, vm_by_vm=vm_by_vm, num_sampling_grid_points=kstar_order, - extra={"velocity_operator": gv_op}, ) # ------------------------------------------------------------------ diff --git a/phono3py/conductivity/protocols.py b/phono3py/conductivity/protocols.py index 968dfd00..0f0376c9 100644 --- a/phono3py/conductivity/protocols.py +++ b/phono3py/conductivity/protocols.py @@ -61,16 +61,23 @@ class VelocityProvider(Protocol): VelocityMatrixProvider Off-diagonal velocity matrix and its outer product (Kubo). - """ + Class attributes + ---------------- + produces_gv_by_gv : bool + True when compute() sets VelocityResult.gv_by_gv. + produces_vm_by_vm : bool + True when compute() sets VelocityResult.vm_by_vm. - def compute(self, gp: GridPointInput) -> VelocityResult: - """Compute velocity quantities at a grid point. + These flags tell the calculator which arrays to pre-allocate + before the grid-point loop. + + """ - The returned VelocityResult must have at minimum - ``group_velocities``, ``gv_by_gv``, and - ``num_sampling_grid_points`` set. + produces_gv_by_gv: bool + produces_vm_by_vm: bool - """ + def compute(self, gp: GridPointInput) -> VelocityResult: + """Compute velocity quantities at a grid point.""" ... @@ -84,19 +91,24 @@ class HeatCapacityProvider(Protocol): HeatCapacityMatrixProvider Heat-capacity matrix Cv_mat (Kubo). + Class attributes + ---------------- + produces_heat_capacity_matrix : bool + True when compute() sets HeatCapacityResult.heat_capacity_matrix. + + These flags tell the calculator which arrays to pre-allocate + before the grid-point loop. + """ + produces_heat_capacity_matrix: bool + def compute( self, gp: GridPointInput, temperatures: NDArray[np.double], ) -> HeatCapacityResult: - """Compute heat-capacity quantities at a grid point. - - The returned HeatCapacityResult must have at minimum - ``heat_capacities`` set (and optionally ``heat_capacity_matrix``). - - """ + """Compute heat-capacity quantities at a grid point.""" ... diff --git a/phono3py/conductivity/rta_calculator.py b/phono3py/conductivity/rta_calculator.py index 20a3df5b..68bad405 100644 --- a/phono3py/conductivity/rta_calculator.py +++ b/phono3py/conductivity/rta_calculator.py @@ -140,7 +140,6 @@ def __init__( self._gamma_U: NDArray[np.double] | None = None self._gamma_elph: NDArray[np.double] | None = None self._gamma_boundary: NDArray[np.double] | None = None - self._extra: dict[str, Any] = {} self._num_sampling_grid_points = 0 self._num_ignored_phonon_modes: NDArray[np.int64] | None = None self._gamma_detail_at_q: NDArray[np.double] | None = None @@ -202,7 +201,6 @@ def run( gamma_elph=self._gamma_elph, vm_by_vm=self._vm_by_vm, heat_capacity_matrix=self._heat_capacity_matrix, - extra=self._extra, ) self._count_ignored_modes(aggregates) self._accumulator.finalize(aggregates) @@ -335,14 +333,14 @@ def get_extra_grid_point_output(self) -> dict[str, Any] | None: """Return per-grid-point extra arrays stored by the calculator. Called by output writers to obtain plugin-defined per-grid-point - arrays (e.g. Wigner velocity_operator) that are written to - the hdf5 file via ``write_kappa_to_hdf5(extra_datasets=...)``. + arrays that are written to the hdf5 file via + ``write_kappa_to_hdf5(extra_datasets=...)``. The caller is responsible for slicing by grid-point index. Returns None when no extra data has been stored. """ - return self._extra if self._extra else None + return None def log_kappa(self) -> None: """Delegate kappa logging to the accumulator. @@ -455,12 +453,16 @@ def _allocate_values(self) -> None: num_temp = len(self._context.temperatures) num_gp = len(self._context.grid_points) num_band0 = len(self._context.band_indices) + num_band = self._context.frequencies.shape[1] self._gamma = np.zeros( (num_sigma, num_temp, num_gp, num_band0), order="C", dtype="double" ) self._gv = np.zeros((num_gp, num_band0, 3), order="C", dtype="double") - self._gv_by_gv = np.zeros((num_gp, num_band0, 6), order="C", dtype="double") + if self._velocity_provider.produces_gv_by_gv: + self._gv_by_gv = np.zeros( + (num_gp, num_band0, 6), order="C", dtype="double" + ) self._cv = np.zeros((num_temp, num_gp, num_band0), order="C", dtype="double") if not self._read_gamma: if self._is_N_U or self._is_gamma_detail: @@ -471,6 +473,22 @@ def _allocate_values(self) -> None: self._gamma_boundary = np.zeros( (num_gp, num_band0), order="C", dtype="double" ) + if self._is_isotope: + self._gamma_iso = np.zeros( + (num_sigma, num_gp, num_band0), order="C", dtype="double" + ) + if self._scattering_provider.is_full_pp: + self._averaged_pp_interaction = np.zeros( + (num_gp, num_band0), order="C", dtype="double" + ) + if self._velocity_provider.produces_vm_by_vm: + self._vm_by_vm = np.zeros( + (num_gp, num_band0, num_band, 6), order="C", dtype="complex128" + ) + if self._cv_provider.produces_heat_capacity_matrix: + self._heat_capacity_matrix = np.zeros( + (num_temp, num_gp, num_band0, num_band), order="C", dtype="double" + ) self._num_ignored_phonon_modes = np.zeros( (num_sigma, num_temp), order="C", dtype="int64" ) @@ -571,19 +589,16 @@ def _run_at_grid_point(self, i_gp: int) -> None: print(" Gamma is read from file.") # Isotope scattering. - gamma_iso: NDArray[np.double] | None = None if self._is_isotope and not self._read_gamma_iso: assert self._isotope_provider is not None gamma_iso = self._isotope_provider.compute(gp_input) - gamma_iso = gamma_iso[:, self._context.band_indices] + self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._context.band_indices] # Boundary scattering (needs group velocities computed above). - gamma_boundary: NDArray[np.double] | None = None if self._boundary_provider is not None: - gamma_boundary = self._boundary_provider.compute( + self._gamma_boundary[i_gp] = self._boundary_provider.compute( vel_result.group_velocities ) - self._gamma_boundary[i_gp] = gamma_boundary if self._log_level: self._show_log(i_gp, vel_result.group_velocities, ave_pp) @@ -591,55 +606,15 @@ def _run_at_grid_point(self, i_gp: int) -> None: # Store per-grid-point output data in calculator. self._gamma[:, :, i_gp, :] = gamma self._gv[i_gp] = vel_result.group_velocities - self._gv_by_gv[i_gp] = vel_result.gv_by_gv + if vel_result.gv_by_gv is not None: + self._gv_by_gv[i_gp] = vel_result.gv_by_gv self._cv[:, i_gp, :] = cv_result.heat_capacities if vel_result.vm_by_vm is not None: - if self._vm_by_vm is None: - self._vm_by_vm = np.zeros( - (len(self._context.grid_points),) + vel_result.vm_by_vm.shape, - dtype="complex128", - order="C", - ) self._vm_by_vm[i_gp] = vel_result.vm_by_vm if cv_result.heat_capacity_matrix is not None: - if self._heat_capacity_matrix is None: - shape = cv_result.heat_capacity_matrix.shape - self._heat_capacity_matrix = np.zeros( - (shape[0], len(self._context.grid_points)) + shape[1:], - dtype="double", - order="C", - ) self._heat_capacity_matrix[:, i_gp] = cv_result.heat_capacity_matrix if ave_pp is not None: - if self._averaged_pp_interaction is None: - self._averaged_pp_interaction = np.zeros( - (len(self._context.grid_points), len(self._context.band_indices)), - dtype="double", - order="C", - ) self._averaged_pp_interaction[i_gp] = ave_pp - if gamma_iso is not None: - if self._gamma_iso is None: - self._gamma_iso = np.zeros( - ( - gamma_iso.shape[0], - len(self._context.grid_points), - gamma_iso.shape[-1], - ), - dtype="double", - order="C", - ) - self._gamma_iso[:, i_gp, :] = gamma_iso - - # Store velocity extra data (e.g. velocity_operator for Wigner). - if vel_result.extra: - for key, val in vel_result.extra.items(): - if key not in self._extra: - num_gp = len(self._context.grid_points) - self._extra[key] = np.zeros( - (num_gp,) + val.shape, dtype=val.dtype, order="C" - ) - self._extra[key][i_gp] = val def _show_log( self, diff --git a/phono3py/conductivity/rta_output.py b/phono3py/conductivity/rta_output.py index f8a7a25a..60463c14 100644 --- a/phono3py/conductivity/rta_output.py +++ b/phono3py/conductivity/rta_output.py @@ -35,7 +35,7 @@ def write_gamma( """Write mode kappa related properties into a hdf5 file.""" grid_points = br.grid_points group_velocities_i = br.group_velocities[i] - gv_by_gv_i = br.gv_by_gv[i] + gv_by_gv_i = br.gv_by_gv[i] if br.gv_by_gv is not None else None extra_gp_full = br.get_extra_grid_point_output() extra_gp_data: dict[str, Any] | None = ( {k: v[i] for k, v in extra_gp_full.items()} diff --git a/phono3py/conductivity/velocity_providers.py b/phono3py/conductivity/velocity_providers.py index feb4bcf5..8ee10ac8 100644 --- a/phono3py/conductivity/velocity_providers.py +++ b/phono3py/conductivity/velocity_providers.py @@ -102,6 +102,9 @@ class GroupVelocityProvider: Verbosity level. Default 0. """ + produces_gv_by_gv: bool = True + produces_vm_by_vm: bool = False + def __init__( self, pp: Interaction, @@ -254,6 +257,9 @@ class VelocityMatrixProvider: """ + produces_gv_by_gv: bool = False + produces_vm_by_vm: bool = True + def __init__( self, pp: Interaction, diff --git a/test/conductivity/ms_smm19/test_phono3py_load_script_ms_smm19.py b/test/conductivity/ms_smm19/test_phono3py_load_script_ms_smm19.py index 9e6432c5..d3c74ff1 100644 --- a/test/conductivity/ms_smm19/test_phono3py_load_script_ms_smm19.py +++ b/test/conductivity/ms_smm19/test_phono3py_load_script_ms_smm19.py @@ -77,11 +77,6 @@ def test_phono3py_load_generates_wigner_kappa_hdf5_contents(): assert np.all(np.isfinite(f["mode_kappa"][:])) assert np.all(np.isfinite(f["mode_kappa_P_RTA"][:])) assert np.all(np.isfinite(f["mode_kappa_C"][:])) - # velocity_operator: (num_gp, num_band0, nat3, 3) - assert "velocity_operator" in f - assert f["velocity_operator"].ndim == 4 - assert f["velocity_operator"].shape[-1] == 3 - assert np.all(np.isfinite(f["velocity_operator"][:])) for created_filename in ( "phono3py.yaml", @@ -221,11 +216,6 @@ def test_phono3py_load_wigner_lbte(): assert np.all(np.isfinite(f["mode_kappa_P_exact"][:])) assert np.all(np.isfinite(f["mode_kappa_P_RTA"][:])) assert np.all(np.isfinite(f["mode_kappa_C"][:])) - # velocity_operator: (num_gp, num_band0, nat3, 3) - assert "velocity_operator" in f - assert f["velocity_operator"].ndim == 4 - assert f["velocity_operator"].shape[-1] == 3 - assert np.all(np.isfinite(f["velocity_operator"][:])) assert kappa_p_exact[0, 0] == pytest.approx(96.14464215062665, abs=0.8) assert kappa_p_rta[0, 0] == pytest.approx(97.30006583719721, abs=0.8) assert kappa_c[0, 0] == pytest.approx(0.1731832844767842, rel=2e-2) @@ -306,17 +296,6 @@ def test_phono3py_load_wigner_write_gamma_contains_isotope_and_N_U(): ) assert np.max(np.abs(gamma - (gamma_n + gamma_u))) < 1e-10 - # Per-grid-point files should contain velocity_operator. - gp_files = sorted(pathlib.Path(".").glob("kappa-m999-g*.hdf5")) - assert gp_files, "No per-grid-point kappa files found" - with h5py.File(gp_files[0], "r") as f_gp: - assert "velocity_operator" in f_gp - vel_op = f_gp["velocity_operator"][:] - # velocity_operator shape: (num_band0, nat3, 3) - assert vel_op.ndim == 3 - assert vel_op.shape[-1] == 3 - assert np.all(np.isfinite(vel_op)) - for created_filename in ( "phono3py.yaml", "fc2.hdf5", From 31929939935a9387899bf5e76ed2ba72f22b70bc Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Tue, 7 Apr 2026 13:23:54 +0900 Subject: [PATCH 2/4] Refactoring of conductivity module --- phono3py/conductivity/grid_point_data.py | 8 +- .../conductivity/heat_capacity_providers.py | 229 +++++++++++++----- phono3py/conductivity/lbte_calculator.py | 90 ++++--- phono3py/conductivity/protocols.py | 17 +- phono3py/conductivity/rta_calculator.py | 157 +++++++----- phono3py/conductivity/scattering_providers.py | 53 ++-- 6 files changed, 343 insertions(+), 211 deletions(-) diff --git a/phono3py/conductivity/grid_point_data.py b/phono3py/conductivity/grid_point_data.py index 62de8f5c..d43ada20 100644 --- a/phono3py/conductivity/grid_point_data.py +++ b/phono3py/conductivity/grid_point_data.py @@ -58,14 +58,14 @@ class VelocityResult: @dataclass class HeatCapacityResult: - """Result from a heat capacity provider at a single grid point. + """Result from a heat capacity provider (bulk computation). Parameters ---------- - heat_capacities : ndarray of double, shape (num_temp, num_band0) - Mode heat capacities (scalar Cv per mode). + heat_capacities : ndarray of double, shape (num_temp, num_gp, num_band0) + Mode heat capacities (scalar Cv per mode) for all grid points. heat_capacity_matrix : ndarray of double, optional - Shape (num_temp, num_band0, num_band). + Shape (num_temp, num_gp, num_band0, num_band). Only set by HeatCapacityMatrixProvider (Kubo). extra : dict Plugin-specific data. diff --git a/phono3py/conductivity/heat_capacity_providers.py b/phono3py/conductivity/heat_capacity_providers.py index e6011bc3..82057b8f 100644 --- a/phono3py/conductivity/heat_capacity_providers.py +++ b/phono3py/conductivity/heat_capacity_providers.py @@ -7,58 +7,145 @@ from phonopy.phonon.thermal_properties import mode_cv from phonopy.physical_units import get_physical_units -from phono3py.conductivity.grid_point_data import GridPointInput, HeatCapacityResult +from phono3py.conductivity.grid_point_data import HeatCapacityResult from phono3py.phonon.heat_capacity_matrix import mode_cv_matrix from phono3py.phonon3.interaction import Interaction +# --------------------------------------------------------------------------- +# Bulk computation functions +# --------------------------------------------------------------------------- -def _get_heat_capacities( - grid_point: int, - pp: Interaction, + +def compute_bulk_mode_cv( + frequencies: NDArray[np.double], + grid_points: NDArray[np.int64], temperatures: NDArray[np.double], + band_indices: NDArray[np.int64], + cutoff_frequency: float, ) -> NDArray[np.double]: - """Return mode heat capacity. - - cv returned should be given to self._cv by + """Compute mode heat capacities for all grid points at once. - self._cv[:, i_data, :] = cv + Parameters + ---------- + frequencies : ndarray of double, shape (num_bz_gp, num_band) + Phonon frequencies in THz for all BZ grid points. + grid_points : ndarray of int64, shape (num_gp,) + BZ grid point indices. + temperatures : ndarray of double, shape (num_temp,) + Temperatures in Kelvin. + band_indices : ndarray of int64, shape (num_band0,) + Selected band indices. + cutoff_frequency : float + Cutoff frequency in THz. + + Returns + ------- + ndarray of double, shape (num_temp, num_gp, num_band0) + Mode heat capacities. """ - if not pp.phonon_all_done: - raise RuntimeError( - "Phonon calculation has not been done yet. " - "Run phono3py.run_phonon_solver() before this method." - ) - - frequencies, _, _ = pp.get_phonons() - assert frequencies is not None - freqs_ev = frequencies[grid_point] * get_physical_units().THzToEv - cv = np.zeros((len(temperatures), len(freqs_ev)), dtype="double") - cutoff = pp.cutoff_frequency * get_physical_units().THzToEv - condition = freqs_ev > cutoff + thz_to_ev = get_physical_units().THzToEv + freqs_ev = frequencies[grid_points] * thz_to_ev # (num_gp, num_band) + cutoff_ev = cutoff_frequency * thz_to_ev + num_gp, num_band = freqs_ev.shape + num_temp = len(temperatures) + flat_freqs = freqs_ev.ravel() # (num_gp * num_band,) with np.errstate(divide="ignore", over="ignore", invalid="ignore"): - cv_vals = mode_cv(temperatures, freqs_ev) - cv = np.where(np.isfinite(cv_vals) & condition[None, :], cv_vals, 0.0) - return cv + cv_flat = mode_cv(temperatures, flat_freqs) # (num_temp, num_gp * num_band) + cv = cv_flat.reshape(num_temp, num_gp, num_band) + + condition = freqs_ev > cutoff_ev # (num_gp, num_band) + cv = np.where(np.isfinite(cv) & condition[None, :, :], cv, 0.0) + return cv[:, :, band_indices] + + +def compute_bulk_cv_matrix( + frequencies: NDArray[np.double], + grid_points: NDArray[np.int64], + temperatures: NDArray[np.double], + band_indices: NDArray[np.int64], + cutoff_frequency: float, +) -> tuple[NDArray[np.double], NDArray[np.double]]: + """Compute cv and cv_matrix for all grid points. + + Parameters + ---------- + frequencies : ndarray of double, shape (num_bz_gp, num_band) + Phonon frequencies in THz for all BZ grid points. + grid_points : ndarray of int64, shape (num_gp,) + BZ grid point indices. + temperatures : ndarray of double, shape (num_temp,) + Temperatures in Kelvin. + band_indices : ndarray of int64, shape (num_band0,) + Selected band indices. + cutoff_frequency : float + Cutoff frequency in THz. + + Returns + ------- + cv : ndarray of double, shape (num_temp, num_gp, num_band0) + Mode heat capacities (diagonal). + cv_mat : ndarray of double, shape (num_temp, num_gp, num_band0, num_band) + Heat capacity matrix. + + """ + thz_to_ev = get_physical_units().THzToEv + num_gp = len(grid_points) + num_band = frequencies.shape[1] + num_temp = len(temperatures) + + cv = compute_bulk_mode_cv( + frequencies, grid_points, temperatures, band_indices, cutoff_frequency + ) + + cutoff_ev = cutoff_frequency * thz_to_ev + cv_mat = np.zeros( + (num_temp, num_gp, len(band_indices), num_band), dtype="double", order="C" + ) + for i_gp, gp in enumerate(grid_points): + freqs_ev = frequencies[gp] * thz_to_ev + condition = freqs_ev > cutoff_ev + condition_2d = np.logical_and.outer(condition, condition) + # cv at this GP for fallback on diagonal + cv_at_gp = compute_bulk_mode_cv( + frequencies, + np.array([gp], dtype="int64"), + temperatures, + np.arange(num_band, dtype="int64"), + cutoff_frequency, + )[:, 0, :] # (num_temp, num_band) + with np.errstate(divide="ignore", over="ignore", invalid="ignore"): + cvm = mode_cv_matrix( + temperatures, freqs_ev + ) # (num_temp, num_band, num_band) + cvm = np.where( + condition_2d[None, :, :], + np.where(np.isfinite(cvm), cvm, cv_at_gp[:, None, :]), + 0.0, + ) + cv_mat[:, i_gp, :, :] = cvm[:, band_indices, :] + + return cv, cv_mat class ModeHeatCapacityProvider: - """Compute scalar mode heat capacities at a grid point. + """Compute scalar mode heat capacities for all grid points at once. This provider implements the ``HeatCapacityProvider`` protocol. It computes the mode heat capacity Cv (per mode, per unit cell) at the requested temperatures using the standard Einstein/harmonic-oscillator - formula, delegating to ``get_heat_capacities`` from ``conductivity.base``. + formula via ``compute_bulk_mode_cv``. The returned ``HeatCapacityResult`` contains ``heat_capacities`` - with shape ``(num_temp, num_band0)``. + with shape ``(num_temp, num_gp, num_band0)``. Parameters ---------- pp : Interaction Interaction instance. Phonon solver must have been run - (``phonon_all_done == True``) before calling ``compute``. + (``phonon_all_done == True``) before calling ``compute_all``. + """ produces_heat_capacity_matrix: bool = False @@ -67,31 +154,43 @@ def __init__(self, pp: Interaction): """Init method.""" self._pp = pp - def compute( + def compute_all( self, - gp: GridPointInput, + frequencies: NDArray[np.double], + grid_points: NDArray[np.int64], temperatures: NDArray[np.double], + band_indices: NDArray[np.int64], + cutoff_frequency: float, ) -> HeatCapacityResult: - """Compute mode heat capacities at a grid point. + """Compute mode heat capacities for all grid points. Parameters ---------- - gp : GridPointInput - Per-grid-point phonon data. + frequencies : ndarray of double, shape (num_bz_gp, num_band) + Phonon frequencies in THz. + grid_points : ndarray of int64, shape (num_gp,) + BZ grid point indices. temperatures : ndarray of double, shape (num_temp,) Temperatures in Kelvin. + band_indices : ndarray of int64, shape (num_band0,) + Selected band indices. + cutoff_frequency : float + Cutoff frequency in THz. Returns ------- HeatCapacityResult - ``heat_capacities`` (num_temp, num_band0) is set. + ``heat_capacities`` (num_temp, num_gp, num_band0) is set. + """ - cv = _get_heat_capacities(gp.grid_point, self._pp, temperatures) - return HeatCapacityResult(heat_capacities=cv[:, self._pp.band_indices]) + cv = compute_bulk_mode_cv( + frequencies, grid_points, temperatures, band_indices, cutoff_frequency + ) + return HeatCapacityResult(heat_capacities=cv) class HeatCapacityMatrixProvider: - """Compute heat capacity matrix at a grid point for the Kubo formula. + """Compute heat capacity matrix for all grid points (Kubo formula). This provider implements the ``HeatCapacityProvider`` protocol for the Green-Kubo formula. It computes the off-diagonal heat capacity matrix @@ -99,16 +198,16 @@ class HeatCapacityMatrixProvider: ``phono3py.phonon.heat_capacity_matrix``. The returned ``HeatCapacityResult`` contains: - - ``heat_capacities`` (num_temp, num_band0): diagonal (standard) mode heat - capacities; set for compatibility with ``RTACalculator``. - - ``heat_capacity_matrix`` (num_temp, num_band0, num_band): full heat - capacity matrix for selected bands (rows) vs all bands (columns). + - ``heat_capacities`` (num_temp, num_gp, num_band0): diagonal (standard) + mode heat capacities. + - ``heat_capacity_matrix`` (num_temp, num_gp, num_band0, num_band): full + heat capacity matrix for selected bands (rows) vs all bands (columns). Parameters ---------- pp : Interaction Interaction instance. Phonon solver must have been run before calling - ``compute``. + ``compute_all``. """ @@ -118,45 +217,41 @@ def __init__(self, pp: Interaction): """Init method.""" self._pp = pp - def compute( + def compute_all( self, - gp: GridPointInput, + frequencies: NDArray[np.double], + grid_points: NDArray[np.int64], temperatures: NDArray[np.double], + band_indices: NDArray[np.int64], + cutoff_frequency: float, ) -> HeatCapacityResult: - """Compute heat capacity matrix at a grid point. + """Compute heat capacity matrix for all grid points. Parameters ---------- - gp : GridPointInput - Per-grid-point phonon data. + frequencies : ndarray of double, shape (num_bz_gp, num_band) + Phonon frequencies in THz. + grid_points : ndarray of int64, shape (num_gp,) + BZ grid point indices. temperatures : ndarray of double, shape (num_temp,) Temperatures in Kelvin. + band_indices : ndarray of int64, shape (num_band0,) + Selected band indices. + cutoff_frequency : float + Cutoff frequency in THz. Returns ------- HeatCapacityResult - ``heat_capacities`` (num_temp, num_band0) and - ``heat_capacity_matrix`` (num_temp, num_band0, num_band) are set. + ``heat_capacities`` (num_temp, num_gp, num_band0) and + ``heat_capacity_matrix`` (num_temp, num_gp, num_band0, num_band) + are set. """ - frequencies, _, _ = self._pp.get_phonons() - assert frequencies is not None - freqs_ev = frequencies[gp.grid_point] * get_physical_units().THzToEv - cv = _get_heat_capacities(gp.grid_point, self._pp, temperatures) - - cutoff = self._pp.cutoff_frequency * get_physical_units().THzToEv - condition = freqs_ev > cutoff - condition = np.logical_and.outer(condition, condition) - - with np.errstate(divide="ignore", over="ignore", invalid="ignore"): - cvm = mode_cv_matrix(temperatures, freqs_ev) - cv_mat = np.where( - condition[None, :, :], - np.where(np.isfinite(cvm), cvm, cv[:, None, :]), - 0.0, - ) - + cv, cv_mat = compute_bulk_cv_matrix( + frequencies, grid_points, temperatures, band_indices, cutoff_frequency + ) return HeatCapacityResult( - heat_capacities=cv[:, self._pp.band_indices], - heat_capacity_matrix=cv_mat[:, self._pp.band_indices, :], + heat_capacities=cv, + heat_capacity_matrix=cv_mat, ) diff --git a/phono3py/conductivity/lbte_calculator.py b/phono3py/conductivity/lbte_calculator.py index f4e138a7..66a9baa1 100644 --- a/phono3py/conductivity/lbte_calculator.py +++ b/phono3py/conductivity/lbte_calculator.py @@ -160,11 +160,20 @@ def run( ) -> None: """Run all grid points and compute kappa. + The computation is organized in separate phases: + + 1. Bulk heat capacity (vectorized, single call). + 2. Velocity loop (per-GP, no interaction state dependency). + 3. Isotope loop (per-GP, no interaction state dependency). + 4. Main collision loop (per-GP, requires Interaction.set_grid_point). + 5. Finalize (assemble collision matrix and compute kappa). + Parameters ---------- on_grid_point : callable or None, optional Called with the grid-point loop index after each grid point is - processed. Used for per-grid-point file writes. + processed in the collision loop. Used for per-grid-point file + writes. """ if self._log_level: @@ -175,11 +184,20 @@ def run( self._prepare_isotope_phonons() - self._num_sampling_grid_points = 0 - self._grid_point_count = 0 + # (1) Bulk heat capacity. + self._compute_bulk_heat_capacities() + + # (2) Velocity loop. + self._compute_all_velocities() + + # (3) Isotope loop. + if self._isotope_provider is not None: + self._compute_all_isotope() + # (4) Main collision loop. + self._grid_point_count = 0 for i_gp in range(len(self._context.ir_grid_points)): - self._run_at_grid_point(i_gp) + self._compute_collision_at_grid_point(i_gp) self._grid_point_count = i_gp + 1 if on_grid_point is not None: on_grid_point(i_gp) @@ -190,6 +208,7 @@ def run( "===================" ) + # (5) Finalize. self._accumulator.finalize(self._build_grid_point_aggregates()) def set_kappa_at_sigmas(self) -> None: @@ -536,43 +555,50 @@ def _build_grid_point_aggregates(self) -> GridPointAggregates: heat_capacity_matrix=self._heat_capacity_matrix, ) - def _run_at_grid_point(self, i_gp: int) -> None: - self._show_log_header(i_gp) - gp_input = self._make_grid_point_input(i_gp) + def _compute_bulk_heat_capacities(self) -> None: + """Compute heat capacities for all grid points at once.""" + cv_result = self._cv_provider.compute_all( + self._context.frequencies, + self._context.ir_grid_points, + self._context.temperatures, + self._context.band_indices, + self._context.cutoff_frequency, + ) + self._cv = cv_result.heat_capacities + if cv_result.heat_capacity_matrix is not None: + self._heat_capacity_matrix = cv_result.heat_capacity_matrix - # Group velocities. - vel_result = self._velocity_provider.compute(gp_input) - gv = vel_result.group_velocities - self._num_sampling_grid_points += vel_result.num_sampling_grid_points + def _compute_all_velocities(self) -> None: + """Compute velocities for all grid points.""" + self._num_sampling_grid_points = 0 + for i_gp in range(len(self._context.ir_grid_points)): + gp_input = self._make_grid_point_input(i_gp) + vel_result = self._velocity_provider.compute(gp_input) + self._num_sampling_grid_points += vel_result.num_sampling_grid_points + self._gv[i_gp] = vel_result.group_velocities + if vel_result.vm_by_vm is not None: + self._vm_by_vm[i_gp] = vel_result.vm_by_vm + + def _compute_all_isotope(self) -> None: + """Compute isotope scattering for all grid points.""" + assert self._isotope_provider is not None + for i_gp in range(len(self._context.ir_grid_points)): + gp_input = self._make_grid_point_input(i_gp) + gamma_iso = self._isotope_provider.compute(gp_input) + self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._context.band_indices] - # Mode heat capacities. - cv_result = self._cv_provider.compute(gp_input, self._context.temperatures) + def _compute_collision_at_grid_point(self, i_gp: int) -> None: + """Compute collision matrix row and gamma at a single grid point.""" + self._show_log_header(i_gp) + gp_input = self._make_grid_point_input(i_gp) - # Collision matrix row + gamma (Stage 1). collision_result = self._collision_provider.compute(gp_input) - - # Store per-grid-point data in calculator. self._gamma[:, :, i_gp, :] = collision_result.gamma - self._gv[i_gp] = gv - self._cv[:, i_gp, :] = cv_result.heat_capacities - # Isotope scattering (optional). - if self._isotope_provider is not None: - gamma_iso = self._isotope_provider.compute(gp_input) - self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._context.band_indices] - - # Averaged pp interaction (optional). if self._is_full_pp and collision_result.averaged_pp is not None: self._averaged_pp_interaction[i_gp] = collision_result.averaged_pp - # Store vm_by_vm and heat_capacity_matrix (for Wigner/Kubo). - if vel_result.vm_by_vm is not None: - self._vm_by_vm[i_gp] = vel_result.vm_by_vm - if cv_result.heat_capacity_matrix is not None: - self._heat_capacity_matrix[:, i_gp] = cv_result.heat_capacity_matrix - - # Accumulate collision matrix row. self._accumulator.store(i_gp, collision_result) if self._log_level: - self._show_log(i_gp, gv) + self._show_log(i_gp, self._gv[i_gp]) diff --git a/phono3py/conductivity/protocols.py b/phono3py/conductivity/protocols.py index 0f0376c9..0f82762e 100644 --- a/phono3py/conductivity/protocols.py +++ b/phono3py/conductivity/protocols.py @@ -12,7 +12,7 @@ VelocityProvider Computes velocity-related quantities at a single BZ grid point. HeatCapacityProvider - Computes heat capacity at a single BZ grid point. + Computes heat capacity for all grid points at once (bulk). ScatteringProvider Computes phonon linewidths at a single BZ grid point. @@ -82,7 +82,7 @@ def compute(self, gp: GridPointInput) -> VelocityResult: class HeatCapacityProvider(Protocol): - """Protocol for computing heat capacity at a grid point. + """Protocol for computing heat capacity for all grid points at once. Built-in implementations ------------------------ @@ -94,7 +94,7 @@ class HeatCapacityProvider(Protocol): Class attributes ---------------- produces_heat_capacity_matrix : bool - True when compute() sets HeatCapacityResult.heat_capacity_matrix. + True when compute_all() sets HeatCapacityResult.heat_capacity_matrix. These flags tell the calculator which arrays to pre-allocate before the grid-point loop. @@ -103,12 +103,15 @@ class HeatCapacityProvider(Protocol): produces_heat_capacity_matrix: bool - def compute( + def compute_all( self, - gp: GridPointInput, + frequencies: NDArray[np.double], + grid_points: NDArray[np.int64], temperatures: NDArray[np.double], + band_indices: NDArray[np.int64], + cutoff_frequency: float, ) -> HeatCapacityResult: - """Compute heat-capacity quantities at a grid point.""" + """Compute heat-capacity quantities for all grid points at once.""" ... @@ -123,7 +126,7 @@ class ScatteringProvider(Protocol): Notes ----- Isotope and boundary scattering are separate diagonal-only contributions - handled by IsotopeScatteringProvider and BoundaryScatteringProvider. + handled by IsotopeScatteringProvider and compute_bulk_boundary_scattering. """ diff --git a/phono3py/conductivity/rta_calculator.py b/phono3py/conductivity/rta_calculator.py index 68bad405..37bec003 100644 --- a/phono3py/conductivity/rta_calculator.py +++ b/phono3py/conductivity/rta_calculator.py @@ -18,9 +18,9 @@ from phono3py.conductivity.heat_capacity_providers import ModeHeatCapacityProvider from phono3py.conductivity.kappa_accumulators import RTAKappaAccumulator from phono3py.conductivity.scattering_providers import ( - BoundaryScatteringProvider, IsotopeScatteringProvider, RTAScatteringProvider, + compute_bulk_boundary_scattering, ) from phono3py.conductivity.utils import ( show_grid_point_frequencies_gv, @@ -112,16 +112,11 @@ def __init__( self._grid_point_count = 0 - # Isotope and boundary providers. + # Isotope provider. self._is_isotope = is_isotope or (mass_variances is not None) self._isotope_provider: IsotopeScatteringProvider | None = None if self._is_isotope: self._isotope_provider = self._build_isotope_provider(mass_variances) - self._boundary_provider: BoundaryScatteringProvider | None = ( - BoundaryScatteringProvider(context.boundary_mfp) - if context.boundary_mfp is not None - else None - ) # Read flags (set via property setters when gamma is loaded from file). self._read_gamma = False @@ -157,11 +152,21 @@ def run( ) -> None: """Run all grid points and compute kappa. + The computation is organized in separate phases: + + 1. Bulk heat capacity (vectorized, single call). + 2. Velocity loop (per-GP, no interaction state dependency). + 3. Isotope loop (per-GP, no interaction state dependency). + 4. Main gamma loop (per-GP, requires Interaction.set_grid_point). + 5. Bulk boundary scattering (vectorized, from accumulated gv). + 6. Finalize (aggregate and compute kappa). + Parameters ---------- on_grid_point : callable or None, optional Called with the grid-point count (0-based index) after each grid - point is processed. Used for per-grid-point file writes. + point is processed in the gamma loop. Used for per-grid-point + file writes. """ if self._context.temperatures is None: @@ -175,11 +180,20 @@ def run( self._prepare_isotope_phonons() - self._num_sampling_grid_points = 0 - self._grid_point_count = 0 + # (1) Bulk heat capacity. + self._compute_bulk_heat_capacities() + # (2) Velocity loop. + self._compute_all_velocities() + + # (3) Isotope loop. + if self._is_isotope and not self._read_gamma_iso: + self._compute_all_isotope() + + # (4) Main gamma loop. + self._grid_point_count = 0 for i_gp in range(len(self._context.grid_points)): - self._run_at_grid_point(i_gp) + self._compute_gamma_at_grid_point(i_gp) self._grid_point_count = i_gp + 1 if on_grid_point is not None: on_grid_point(i_gp) @@ -190,18 +204,14 @@ def run( "===================" ) - aggregates = GridPointAggregates( - num_sampling_grid_points=self._num_sampling_grid_points, - group_velocities=self._gv, - mode_heat_capacities=self._cv, - gv_by_gv=self._gv_by_gv, - gamma=self._gamma, - gamma_isotope=self._gamma_iso, - gamma_boundary=self._gamma_boundary, - gamma_elph=self._gamma_elph, - vm_by_vm=self._vm_by_vm, - heat_capacity_matrix=self._heat_capacity_matrix, - ) + # (5) Bulk boundary scattering. + if self._context.boundary_mfp is not None: + self._gamma_boundary = compute_bulk_boundary_scattering( + self._gv, self._context.boundary_mfp + ) + + # (6) Finalize. + aggregates = self._build_grid_point_aggregates() self._count_ignored_modes(aggregates) self._accumulator.finalize(aggregates) @@ -460,19 +470,13 @@ def _allocate_values(self) -> None: ) self._gv = np.zeros((num_gp, num_band0, 3), order="C", dtype="double") if self._velocity_provider.produces_gv_by_gv: - self._gv_by_gv = np.zeros( - (num_gp, num_band0, 6), order="C", dtype="double" - ) + self._gv_by_gv = np.zeros((num_gp, num_band0, 6), order="C", dtype="double") self._cv = np.zeros((num_temp, num_gp, num_band0), order="C", dtype="double") if not self._read_gamma: if self._is_N_U or self._is_gamma_detail: shape = (num_sigma, num_temp, num_gp, num_band0) self._gamma_N = np.zeros(shape, order="C", dtype="double") self._gamma_U = np.zeros(shape, order="C", dtype="double") - if self._context.boundary_mfp is not None: - self._gamma_boundary = np.zeros( - (num_gp, num_band0), order="C", dtype="double" - ) if self._is_isotope: self._gamma_iso = np.zeros( (num_sigma, num_gp, num_band0), order="C", dtype="double" @@ -557,20 +561,47 @@ def _show_log_header(self, i_gp: int) -> None: mass_variances=mass_variances, ) - def _run_at_grid_point(self, i_gp: int) -> None: - self._show_log_header(i_gp) - gp_input = self._make_grid_point_input(i_gp) + def _compute_bulk_heat_capacities(self) -> None: + """Compute heat capacities for all grid points at once.""" + assert self._context.temperatures is not None + cv_result = self._cv_provider.compute_all( + self._context.frequencies, + self._context.grid_points, + self._context.temperatures, + self._context.band_indices, + self._context.cutoff_frequency, + ) + self._cv = cv_result.heat_capacities + if cv_result.heat_capacity_matrix is not None: + self._heat_capacity_matrix = cv_result.heat_capacity_matrix - # Velocity. - vel_result = self._velocity_provider.compute(gp_input) - self._num_sampling_grid_points += vel_result.num_sampling_grid_points + def _compute_all_velocities(self) -> None: + """Compute velocities for all grid points.""" + self._num_sampling_grid_points = 0 + for i_gp in range(len(self._context.grid_points)): + gp_input = self._make_grid_point_input(i_gp) + vel_result = self._velocity_provider.compute(gp_input) + self._num_sampling_grid_points += vel_result.num_sampling_grid_points + self._gv[i_gp] = vel_result.group_velocities + if vel_result.gv_by_gv is not None: + self._gv_by_gv[i_gp] = vel_result.gv_by_gv + if vel_result.vm_by_vm is not None: + self._vm_by_vm[i_gp] = vel_result.vm_by_vm + + def _compute_all_isotope(self) -> None: + """Compute isotope scattering for all grid points.""" + assert self._isotope_provider is not None + for i_gp in range(len(self._context.grid_points)): + gp_input = self._make_grid_point_input(i_gp) + gamma_iso = self._isotope_provider.compute(gp_input) + self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._context.band_indices] - # Heat capacity. - assert self._context.temperatures is not None - cv_result = self._cv_provider.compute(gp_input, self._context.temperatures) + def _compute_gamma_at_grid_point(self, i_gp: int) -> None: + """Compute ph-ph scattering gamma at a single grid point.""" + self._show_log_header(i_gp) - # ph-ph scattering. if not self._read_gamma: + gp_input = self._make_grid_point_input(i_gp) scat_result = self._scattering_provider.compute(gp_input) gamma = scat_result.gamma ave_pp = scat_result.averaged_pp_interaction @@ -582,39 +613,31 @@ def _run_at_grid_point(self, i_gp: int) -> None: if g_U is not None and self._gamma_U is not None: self._gamma_U[:, :, i_gp, :] = g_U self._gamma_detail_at_q = self._scattering_provider.gamma_detail_at_q + self._gamma[:, :, i_gp, :] = gamma + if ave_pp is not None and self._averaged_pp_interaction is not None: + self._averaged_pp_interaction[i_gp] = ave_pp else: - gamma = self._gamma[:, :, i_gp, :] ave_pp = None if self._log_level: print(" Gamma is read from file.") - # Isotope scattering. - if self._is_isotope and not self._read_gamma_iso: - assert self._isotope_provider is not None - gamma_iso = self._isotope_provider.compute(gp_input) - self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._context.band_indices] - - # Boundary scattering (needs group velocities computed above). - if self._boundary_provider is not None: - self._gamma_boundary[i_gp] = self._boundary_provider.compute( - vel_result.group_velocities - ) - if self._log_level: - self._show_log(i_gp, vel_result.group_velocities, ave_pp) - - # Store per-grid-point output data in calculator. - self._gamma[:, :, i_gp, :] = gamma - self._gv[i_gp] = vel_result.group_velocities - if vel_result.gv_by_gv is not None: - self._gv_by_gv[i_gp] = vel_result.gv_by_gv - self._cv[:, i_gp, :] = cv_result.heat_capacities - if vel_result.vm_by_vm is not None: - self._vm_by_vm[i_gp] = vel_result.vm_by_vm - if cv_result.heat_capacity_matrix is not None: - self._heat_capacity_matrix[:, i_gp] = cv_result.heat_capacity_matrix - if ave_pp is not None: - self._averaged_pp_interaction[i_gp] = ave_pp + self._show_log(i_gp, self._gv[i_gp], ave_pp) + + def _build_grid_point_aggregates(self) -> GridPointAggregates: + """Build GridPointAggregates for accumulator.finalize().""" + return GridPointAggregates( + num_sampling_grid_points=self._num_sampling_grid_points, + group_velocities=self._gv, + mode_heat_capacities=self._cv, + gv_by_gv=self._gv_by_gv, + gamma=self._gamma, + gamma_isotope=self._gamma_iso, + gamma_boundary=self._gamma_boundary, + gamma_elph=self._gamma_elph, + vm_by_vm=self._vm_by_vm, + heat_capacity_matrix=self._heat_capacity_matrix, + ) def _show_log( self, diff --git a/phono3py/conductivity/scattering_providers.py b/phono3py/conductivity/scattering_providers.py index e10219a4..257249e8 100644 --- a/phono3py/conductivity/scattering_providers.py +++ b/phono3py/conductivity/scattering_providers.py @@ -471,50 +471,35 @@ def compute(self, gp: GridPointInput) -> NDArray[np.double]: return np.array(gamma_iso, dtype="double", order="C") -class BoundaryScatteringProvider: - """Compute boundary scattering linewidth at a grid point. +def compute_bulk_boundary_scattering( + group_velocities: NDArray[np.double], + boundary_mfp: float, +) -> NDArray[np.double]: + """Compute boundary scattering linewidth for all grid points at once. - Returns gamma_boundary as NDArray with shape ``(num_band0,)``. The formula is: - gamma_boundary[s] = |v_s| * 1e6 * Angstrom / (4 * pi * boundary_mfp) + gamma_boundary[gp, s] = |v_s| * 1e6 * Angstrom / (4 * pi * boundary_mfp) where ``boundary_mfp`` is in micrometres and ``|v_s|`` is the group velocity magnitude in THz*Angstrom. Parameters ---------- + group_velocities : ndarray of double, shape (num_gp, num_band0, 3) + Group velocities in THz*Angstrom. boundary_mfp : float Boundary mean free path in micrometres. - """ - - def __init__(self, boundary_mfp: float): - """Init method.""" - self._boundary_mfp = boundary_mfp - def compute( - self, - group_velocities: NDArray[np.double], - ) -> NDArray[np.double]: - """Compute boundary scattering linewidth. - - Parameters - ---------- - group_velocities : ndarray, shape (num_band0, 3) - Group velocities in THz*Angstrom. + Returns + ------- + ndarray of double, shape (num_gp, num_band0) + Boundary scattering linewidth. - Returns - ------- - ndarray of double, shape (num_band0,) - Boundary scattering linewidth. - """ - num_band0 = len(group_velocities) - g_boundary = np.zeros(num_band0, dtype="double") - for ll in range(num_band0): - g_boundary[ll] = ( - np.linalg.norm(group_velocities[ll]) - * get_physical_units().Angstrom - * 1e6 - / (4 * np.pi * self._boundary_mfp) - ) - return g_boundary + """ + return ( + np.linalg.norm(group_velocities, axis=-1) + * get_physical_units().Angstrom + * 1e6 + / (4 * np.pi * boundary_mfp) + ) From b04f91c6205a4fa0438034a9bf68ea166a12b833 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Tue, 7 Apr 2026 17:28:26 +0900 Subject: [PATCH 3/4] Refactoring of conductivity module --- phono3py/api_phono3py.py | 4 +-- phono3py/conductivity/lbte_calculator.py | 20 +++++++---- phono3py/conductivity/lbte_init.py | 11 +++--- phono3py/conductivity/rta_calculator.py | 34 +++++++++++-------- phono3py/conductivity/scattering_providers.py | 4 --- phono3py/cui/phono3py_script.py | 7 ++++ phono3py/cui/settings.py | 2 +- 7 files changed, 49 insertions(+), 33 deletions(-) diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index be8503bb..e8a59937 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -2248,7 +2248,7 @@ def run_thermal_conductivity( is_kappa_star: bool = True, gv_delta_q: float | None = None, # for group velocity is_full_pp: bool = False, - pinv_cutoff: float = 1.0e-8, # for pseudo-inversion of collision matrix + pinv_cutoff: float | None = None, # for pseudo-inversion of collision matrix pinv_method: int = 0, # for pseudo-inversion of collision matrix pinv_solver: int = 0, # solver of pseudo-inversion of collision matrix write_gamma: bool = False, @@ -2320,7 +2320,7 @@ def run_thermal_conductivity( improve of efficiency is expected. With smearing method, even if this is set False, full elements are computed unless `sigma_cutoff` is specified. - pinv_cutoff : float, optional, default is 1.0e-8 + pinv_cutoff : float, optional, default is None (usually 1.0e-8) Direct solution only (`is_LBTE=True`). This is used as a criterion to judge the eigenvalues are considered as zero or not in pseudo-inversion of collision matrix. See also `pinv_method`. diff --git a/phono3py/conductivity/lbte_calculator.py b/phono3py/conductivity/lbte_calculator.py index 66a9baa1..973411a1 100644 --- a/phono3py/conductivity/lbte_calculator.py +++ b/phono3py/conductivity/lbte_calculator.py @@ -176,22 +176,28 @@ def run( writes. """ - if self._log_level: - print( - "==================== Lattice thermal conductivity (LBTE) " - "====================" - ) - self._prepare_isotope_phonons() # (1) Bulk heat capacity. - self._compute_bulk_heat_capacities() + if self._log_level: + print("Running heat capacity calculations...") + self._compute_bulk_heat_capacities() # (2) Velocity loop. + if self._log_level: + print("Running velocity calculations...") self._compute_all_velocities() # (3) Isotope loop. if self._isotope_provider is not None: + if self._log_level: + for sigma in self._context.sigmas: + print("Running isotope scattering calculations ", end="") + print( + "with tetrahedron method..." + if sigma is None + else f"sigma={sigma}..." + ) self._compute_all_isotope() # (4) Main collision loop. diff --git a/phono3py/conductivity/lbte_init.py b/phono3py/conductivity/lbte_init.py index 699f2324..e5d6e99f 100644 --- a/phono3py/conductivity/lbte_init.py +++ b/phono3py/conductivity/lbte_init.py @@ -173,7 +173,7 @@ def get_thermal_conductivity_LBTE( gv_delta_q: float | None = None, is_full_pp: bool = False, transport_type: str | None = None, - pinv_cutoff: float = 1.0e-8, + pinv_cutoff: float | None = None, pinv_solver: int = 0, # default: dsyev in lapacke pinv_method: int = 0, # default: abs(eig) < cutoff write_collision: bool = False, @@ -197,14 +197,15 @@ def get_thermal_conductivity_LBTE( _grid_points = ( np.asarray(grid_points, dtype="int64") if grid_points is not None else None ) + if pinv_cutoff is None: + _pinv_cutoff = 1.0e-8 + else: + _pinv_cutoff = pinv_cutoff if sigmas is None: sigmas = [] if log_level: print("-" * 19 + " Lattice thermal conductivity (LBTE) " + "-" * 19) - print( - "Cutoff frequency of pseudo inversion of collision matrix: %s" % pinv_cutoff - ) method = f"{transport_type}-lbte" if transport_type else "lbte" return _run_standard_lbte( @@ -222,7 +223,7 @@ def get_thermal_conductivity_LBTE( is_kappa_star=is_kappa_star, gv_delta_q=gv_delta_q, is_full_pp=is_full_pp, - pinv_cutoff=pinv_cutoff, + pinv_cutoff=_pinv_cutoff, pinv_solver=pinv_solver, pinv_method=pinv_method, write_collision=write_collision, diff --git a/phono3py/conductivity/rta_calculator.py b/phono3py/conductivity/rta_calculator.py index 37bec003..8637c6c1 100644 --- a/phono3py/conductivity/rta_calculator.py +++ b/phono3py/conductivity/rta_calculator.py @@ -172,25 +172,37 @@ def run( if self._context.temperatures is None: raise RuntimeError("Set temperatures before calling run().") - if self._log_level: - print( - "==================== Lattice thermal conductivity (RTA) " - "====================" - ) - self._prepare_isotope_phonons() # (1) Bulk heat capacity. + if self._log_level: + print("Running heat capacity calculations...") self._compute_bulk_heat_capacities() # (2) Velocity loop. + if self._log_level: + print("Running velocity calculations...") self._compute_all_velocities() - # (3) Isotope loop. + # (3) Bulk boundary scattering. + if self._context.boundary_mfp is not None: + self._gamma_boundary = compute_bulk_boundary_scattering( + self._gv, self._context.boundary_mfp + ) + + # (4) Isotope loop. if self._is_isotope and not self._read_gamma_iso: + if self._log_level: + for sigma in self._context.sigmas: + print("Running isotope scattering calculations ", end="") + print( + "with tetrahedron method..." + if sigma is None + else f"sigma={sigma}..." + ) self._compute_all_isotope() - # (4) Main gamma loop. + # (5) Main gamma loop. self._grid_point_count = 0 for i_gp in range(len(self._context.grid_points)): self._compute_gamma_at_grid_point(i_gp) @@ -204,12 +216,6 @@ def run( "===================" ) - # (5) Bulk boundary scattering. - if self._context.boundary_mfp is not None: - self._gamma_boundary = compute_bulk_boundary_scattering( - self._gv, self._context.boundary_mfp - ) - # (6) Finalize. aggregates = self._build_grid_point_aggregates() self._count_ignored_modes(aggregates) diff --git a/phono3py/conductivity/scattering_providers.py b/phono3py/conductivity/scattering_providers.py index 257249e8..328afd7f 100644 --- a/phono3py/conductivity/scattering_providers.py +++ b/phono3py/conductivity/scattering_providers.py @@ -459,10 +459,6 @@ def compute(self, gp: GridPointInput) -> NDArray[np.double]: """ gamma_iso = [] for sigma in self._sigmas: - if self._log_level: - text = "Calculating Gamma of ph-isotope with " - text += "tetrahedron method" if sigma is None else "sigma=%s" % sigma - print(text) self._isotope.sigma = sigma self._isotope.set_grid_point(gp.grid_point) self._isotope.run() diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py index 7ad6df0b..bc5f61f8 100644 --- a/phono3py/cui/phono3py_script.py +++ b/phono3py/cui/phono3py_script.py @@ -1027,6 +1027,13 @@ def _init_phph_interaction( else: print("fc3-r2q-transformation over three atoms: False") + if settings.transport_type is not None: + print(f"Optional transport type: {settings.transport_type.upper()}") + if settings.pinv_cutoff is not None: + print( + f"Pseudo inversion cutoff of collision matrix: {settings.pinv_cutoff}" + ) + ave_pp = settings.constant_averaged_pp_interaction phono3py.init_phph_interaction( nac_q_direction=settings.nac_q_direction, diff --git a/phono3py/cui/settings.py b/phono3py/cui/settings.py index 8168fd90..cfffafaf 100644 --- a/phono3py/cui/settings.py +++ b/phono3py/cui/settings.py @@ -96,7 +96,7 @@ def __init__(self, load_phono3py_yaml: bool = False) -> None: self.read_pp: bool = False self.output_yaml_filename: str | os.PathLike | None = None self.phonon_supercell_matrix: NDArray[np.int64] | None = None - self.pinv_cutoff: float = 1.0e-8 + self.pinv_cutoff: float | None = None self.pinv_solver: int = 0 self.pinv_method: int = 0 self.pp_conversion_factor: float | None = None From 70ff05f46bbcfaf74a125ac71f8850d89957d4aa Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Tue, 7 Apr 2026 19:01:39 +0900 Subject: [PATCH 4/4] Refactoring of conductivity module --- phono3py/conductivity/__init__.py | 3 - phono3py/conductivity/build_components.py | 4 +- phono3py/conductivity/factory.py | 12 ++- phono3py/conductivity/grid_point_data.py | 94 +------------------ .../conductivity/heat_capacity_providers.py | 52 ++++------ phono3py/conductivity/lbte_calculator.py | 36 ++----- .../conductivity/lbte_collision_provider.py | 27 ++---- .../ms_smm19/velocity_providers.py | 34 +++---- phono3py/conductivity/njc23/__init__.py | 4 +- phono3py/conductivity/protocols.py | 19 ++-- phono3py/conductivity/rta_calculator.py | 32 ++----- phono3py/conductivity/scattering_providers.py | 45 +++++---- phono3py/conductivity/smm19/__init__.py | 4 +- phono3py/conductivity/velocity_providers.py | 59 ++++++------ phono3py/phonon/grid.py | 2 +- test/conductivity/test_grid_point_data.py | 65 +++++-------- 16 files changed, 161 insertions(+), 331 deletions(-) diff --git a/phono3py/conductivity/__init__.py b/phono3py/conductivity/__init__.py index 54dc7e4a..fd71a60e 100644 --- a/phono3py/conductivity/__init__.py +++ b/phono3py/conductivity/__init__.py @@ -6,7 +6,6 @@ register_variant, VariantBuildContext, GridPointAggregates, - GridPointInput, VelocityResult, HeatCapacityResult, ScatteringResult, @@ -20,7 +19,6 @@ from phono3py.conductivity.factory import register_variant from phono3py.conductivity.grid_point_data import ( GridPointAggregates, - GridPointInput, HeatCapacityResult, ScatteringResult, VelocityResult, @@ -35,7 +33,6 @@ "register_variant", "VariantBuildContext", "GridPointAggregates", - "GridPointInput", "VelocityResult", "HeatCapacityResult", "ScatteringResult", diff --git a/phono3py/conductivity/build_components.py b/phono3py/conductivity/build_components.py index 2a2cf700..d0ff01b0 100644 --- a/phono3py/conductivity/build_components.py +++ b/phono3py/conductivity/build_components.py @@ -417,10 +417,10 @@ def _resolve_rta_grid_points( return np.array(grid_points, dtype="int64"), ir_gps_bzg, gp_weights if not is_kappa_star: - all_gps = np.array(bz_grid.grg2bzg, dtype="int64") + all_gps = bz_grid.grg2bzg return all_gps, all_gps, np.ones(len(all_gps), dtype="int64") - return ir_gps_bzg, ir_gps_bzg, np.array(ir_weights, dtype="int64") + return ir_gps_bzg, ir_gps_bzg, ir_weights def build_rta_base_components( diff --git a/phono3py/conductivity/factory.py b/phono3py/conductivity/factory.py index 3cf35fb0..4e7bb07e 100644 --- a/phono3py/conductivity/factory.py +++ b/phono3py/conductivity/factory.py @@ -143,7 +143,7 @@ def _rta_factory( cv_provider = ( make_cv_provider(ctx) if make_cv_provider is not None - else ModeHeatCapacityProvider(interaction) + else ModeHeatCapacityProvider(interaction, base.context.temperatures) ) accumulator = make_rta_accumulator(ctx) @@ -189,7 +189,7 @@ def _lbte_factory( cv_provider = ( make_cv_provider(ctx) if make_cv_provider is not None - else ModeHeatCapacityProvider(interaction) + else ModeHeatCapacityProvider(interaction, base.context.temperatures) ) accumulator = make_lbte_accumulator(ctx) @@ -371,11 +371,13 @@ def _make_rta_calculator( interaction, point_operations=base.point_ops, rotations_cartesian=base.rot_cart, + grid_points=base.context.grid_points, + grid_weights=base.context.grid_weights, is_kappa_star=config.is_kappa_star, gv_delta_q=config.gv_delta_q, log_level=config.log_level, ) - cv_provider = ModeHeatCapacityProvider(interaction) + cv_provider = ModeHeatCapacityProvider(interaction, base.context.temperatures) accumulator = RTAKappaAccumulator( context=base.context, conversion_factor=conversion_factor, @@ -406,11 +408,13 @@ def _make_lbte_calculator( interaction, point_operations=base.point_ops, rotations_cartesian=base.rot_cart, + grid_points=base.context.ir_grid_points, + grid_weights=base.context.grid_weights, is_kappa_star=config.is_kappa_star, gv_delta_q=config.gv_delta_q, log_level=config.log_level, ) - cv_provider = ModeHeatCapacityProvider(interaction) + cv_provider = ModeHeatCapacityProvider(interaction, base.context.temperatures) accumulator = LBTEKappaAccumulator( base.solver, context=base.context, log_level=config.log_level ) diff --git a/phono3py/conductivity/grid_point_data.py b/phono3py/conductivity/grid_point_data.py index d43ada20..4c1ddae1 100644 --- a/phono3py/conductivity/grid_point_data.py +++ b/phono3py/conductivity/grid_point_data.py @@ -1,8 +1,8 @@ """Per-grid-point data containers. -This module defines GridPointInput, provider result types -(VelocityResult, HeatCapacityResult, ScatteringResult), and -GridPointAggregates used in the conductivity calculation. +This module defines provider result types (VelocityResult, +HeatCapacityResult, ScatteringResult) and GridPointAggregates used in +the conductivity calculation. Protocol interfaces (VelocityProvider, HeatCapacityProvider, ScatteringProvider) are defined in ``phono3py.conductivity.protocols`` @@ -17,8 +17,6 @@ import numpy as np from numpy.typing import NDArray -from phono3py.phonon.grid import BZGrid, get_qpoints_from_bz_grid_points - # --------------------------------------------------------------------------- # Per-grid-point data containers # --------------------------------------------------------------------------- @@ -97,42 +95,6 @@ class ScatteringResult: extra: dict[str, Any] = field(default_factory=dict) -# --------------------------------------------------------------------------- -# Grid-point input / result containers -# --------------------------------------------------------------------------- - - -@dataclass -class GridPointInput: - """Phonon quantities at a single irreducible BZ grid point. - - All fields are filled before being passed to building blocks. - - Parameters - ---------- - grid_point : int - BZ grid point index. - q_point : ndarray, shape (3,) - q-point coordinates in reduced reciprocal coordinates. - frequencies : ndarray, shape (num_band,) - Phonon frequencies in THz for *all* bands at this grid point. - eigenvectors : ndarray, shape (num_band, num_band), complex - Phonon eigenvectors at this grid point. - grid_weight : int - Symmetry weight for BZ summation (number of arms of the k-star). - band_indices : ndarray of int64, shape (num_band0,) - Selected band indices. num_band0 <= num_band. - - """ - - grid_point: int - q_point: NDArray[np.double] - frequencies: NDArray[np.double] - eigenvectors: NDArray[np.cdouble] - grid_weight: int - band_indices: NDArray[np.int64] - - # --------------------------------------------------------------------------- # Aggregated per-grid-point data (Calculator -> Accumulator.finalize()) # --------------------------------------------------------------------------- @@ -229,54 +191,6 @@ def compute_effective_gamma( return out -def make_grid_point_input( - grid_point: int, - grid_weight: int, - frequencies: NDArray[np.double], - eigenvectors: NDArray[np.cdouble], - bz_grid: BZGrid, - band_indices: NDArray[np.int64], -) -> GridPointInput: - """Create a GridPointInput for a single BZ grid point. - - Parameters - ---------- - grid_point : int - BZ grid point index. - grid_weight : int - Symmetry weight for BZ summation. - frequencies : ndarray of double, shape (num_bz_gp, num_band) - Phonon frequencies array indexed by BZ grid point. - eigenvectors : ndarray of cdouble, shape (num_bz_gp, num_band, num_band) - Phonon eigenvectors array indexed by BZ grid point. - bz_grid : BZGrid - Brillouin zone grid object. - band_indices : ndarray of int64, shape (num_band0,) - Selected band indices. - - Returns - ------- - GridPointInput - - """ - return GridPointInput( - grid_point=grid_point, - q_point=np.array( - get_qpoints_from_bz_grid_points(grid_point, bz_grid), - dtype="double", - ), - frequencies=frequencies[grid_point], - eigenvectors=eigenvectors[grid_point], - grid_weight=grid_weight, - band_indices=band_indices, - ) - - -# --------------------------------------------------------------------------- -# Building-block Protocol interfaces -# --------------------------------------------------------------------------- - - # --------------------------------------------------------------------------- # Re-export Protocol interfaces for backward compatibility. # Canonical definitions live in phono3py.conductivity.protocols. @@ -290,12 +204,10 @@ def make_grid_point_input( __all__ = [ "GridPointAggregates", - "GridPointInput", "HeatCapacityResult", "ScatteringResult", "VelocityResult", "compute_effective_gamma", - "make_grid_point_input", "VelocityProvider", "HeatCapacityProvider", "ScatteringProvider", diff --git a/phono3py/conductivity/heat_capacity_providers.py b/phono3py/conductivity/heat_capacity_providers.py index 82057b8f..5e0224d2 100644 --- a/phono3py/conductivity/heat_capacity_providers.py +++ b/phono3py/conductivity/heat_capacity_providers.py @@ -144,38 +144,27 @@ class ModeHeatCapacityProvider: ---------- pp : Interaction Interaction instance. Phonon solver must have been run - (``phonon_all_done == True``) before calling ``compute_all``. + (``phonon_all_done == True``) before calling ``compute``. """ produces_heat_capacity_matrix: bool = False - def __init__(self, pp: Interaction): + def __init__(self, pp: Interaction, temperatures: NDArray[np.double]): """Init method.""" self._pp = pp + self._temperatures = temperatures - def compute_all( + def compute( self, - frequencies: NDArray[np.double], grid_points: NDArray[np.int64], - temperatures: NDArray[np.double], - band_indices: NDArray[np.int64], - cutoff_frequency: float, ) -> HeatCapacityResult: """Compute mode heat capacities for all grid points. Parameters ---------- - frequencies : ndarray of double, shape (num_bz_gp, num_band) - Phonon frequencies in THz. grid_points : ndarray of int64, shape (num_gp,) BZ grid point indices. - temperatures : ndarray of double, shape (num_temp,) - Temperatures in Kelvin. - band_indices : ndarray of int64, shape (num_band0,) - Selected band indices. - cutoff_frequency : float - Cutoff frequency in THz. Returns ------- @@ -183,8 +172,13 @@ def compute_all( ``heat_capacities`` (num_temp, num_gp, num_band0) is set. """ + frequencies = self._pp.get_phonons()[0] cv = compute_bulk_mode_cv( - frequencies, grid_points, temperatures, band_indices, cutoff_frequency + frequencies, + grid_points, + self._temperatures, + self._pp.band_indices, + self._pp.cutoff_frequency, ) return HeatCapacityResult(heat_capacities=cv) @@ -207,38 +201,27 @@ class HeatCapacityMatrixProvider: ---------- pp : Interaction Interaction instance. Phonon solver must have been run before calling - ``compute_all``. + ``compute``. """ produces_heat_capacity_matrix: bool = True - def __init__(self, pp: Interaction): + def __init__(self, pp: Interaction, temperatures: NDArray[np.double]): """Init method.""" self._pp = pp + self._temperatures = temperatures - def compute_all( + def compute( self, - frequencies: NDArray[np.double], grid_points: NDArray[np.int64], - temperatures: NDArray[np.double], - band_indices: NDArray[np.int64], - cutoff_frequency: float, ) -> HeatCapacityResult: """Compute heat capacity matrix for all grid points. Parameters ---------- - frequencies : ndarray of double, shape (num_bz_gp, num_band) - Phonon frequencies in THz. grid_points : ndarray of int64, shape (num_gp,) BZ grid point indices. - temperatures : ndarray of double, shape (num_temp,) - Temperatures in Kelvin. - band_indices : ndarray of int64, shape (num_band0,) - Selected band indices. - cutoff_frequency : float - Cutoff frequency in THz. Returns ------- @@ -248,8 +231,13 @@ def compute_all( are set. """ + frequencies = self._pp.get_phonons()[0] cv, cv_mat = compute_bulk_cv_matrix( - frequencies, grid_points, temperatures, band_indices, cutoff_frequency + frequencies, + grid_points, + self._temperatures, + self._pp.band_indices, + self._pp.cutoff_frequency, ) return HeatCapacityResult( heat_capacities=cv, diff --git a/phono3py/conductivity/lbte_calculator.py b/phono3py/conductivity/lbte_calculator.py index 973411a1..86a64ab3 100644 --- a/phono3py/conductivity/lbte_calculator.py +++ b/phono3py/conductivity/lbte_calculator.py @@ -43,11 +43,7 @@ from numpy.typing import NDArray from phono3py.conductivity.context import ConductivityContext -from phono3py.conductivity.grid_point_data import ( - GridPointAggregates, - GridPointInput, - make_grid_point_input, -) +from phono3py.conductivity.grid_point_data import GridPointAggregates from phono3py.conductivity.heat_capacity_providers import ModeHeatCapacityProvider from phono3py.conductivity.kappa_accumulators import LBTEKappaAccumulator from phono3py.conductivity.lbte_collision_provider import LBTECollisionProvider @@ -481,16 +477,6 @@ def _prepare_isotope_phonons(self) -> None: # Private: per-grid-point computation # ------------------------------------------------------------------ - def _make_grid_point_input(self, i_gp: int) -> GridPointInput: - return make_grid_point_input( - grid_point=int(self._context.ir_grid_points[i_gp]), - grid_weight=int(self._context.grid_weights[i_gp]), - frequencies=self._context.frequencies, - eigenvectors=self._context.eigenvectors, - bz_grid=self._context.bz_grid, - band_indices=self._context.band_indices, - ) - def _show_log_header(self, i_gp: int) -> None: if not self._log_level: return @@ -563,13 +549,7 @@ def _build_grid_point_aggregates(self) -> GridPointAggregates: def _compute_bulk_heat_capacities(self) -> None: """Compute heat capacities for all grid points at once.""" - cv_result = self._cv_provider.compute_all( - self._context.frequencies, - self._context.ir_grid_points, - self._context.temperatures, - self._context.band_indices, - self._context.cutoff_frequency, - ) + cv_result = self._cv_provider.compute(self._context.ir_grid_points) self._cv = cv_result.heat_capacities if cv_result.heat_capacity_matrix is not None: self._heat_capacity_matrix = cv_result.heat_capacity_matrix @@ -578,8 +558,8 @@ def _compute_all_velocities(self) -> None: """Compute velocities for all grid points.""" self._num_sampling_grid_points = 0 for i_gp in range(len(self._context.ir_grid_points)): - gp_input = self._make_grid_point_input(i_gp) - vel_result = self._velocity_provider.compute(gp_input) + grid_point = int(self._context.ir_grid_points[i_gp]) + vel_result = self._velocity_provider.compute(grid_point) self._num_sampling_grid_points += vel_result.num_sampling_grid_points self._gv[i_gp] = vel_result.group_velocities if vel_result.vm_by_vm is not None: @@ -589,16 +569,16 @@ def _compute_all_isotope(self) -> None: """Compute isotope scattering for all grid points.""" assert self._isotope_provider is not None for i_gp in range(len(self._context.ir_grid_points)): - gp_input = self._make_grid_point_input(i_gp) - gamma_iso = self._isotope_provider.compute(gp_input) + grid_point = int(self._context.ir_grid_points[i_gp]) + gamma_iso = self._isotope_provider.compute(grid_point) self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._context.band_indices] def _compute_collision_at_grid_point(self, i_gp: int) -> None: """Compute collision matrix row and gamma at a single grid point.""" self._show_log_header(i_gp) - gp_input = self._make_grid_point_input(i_gp) + grid_point = int(self._context.ir_grid_points[i_gp]) - collision_result = self._collision_provider.compute(gp_input) + collision_result = self._collision_provider.compute(grid_point) self._gamma[:, :, i_gp, :] = collision_result.gamma if self._is_full_pp and collision_result.averaged_pp is not None: diff --git a/phono3py/conductivity/lbte_collision_provider.py b/phono3py/conductivity/lbte_collision_provider.py index dcc90162..b7c9b507 100644 --- a/phono3py/conductivity/lbte_collision_provider.py +++ b/phono3py/conductivity/lbte_collision_provider.py @@ -42,7 +42,6 @@ import numpy as np from numpy.typing import NDArray -from phono3py.conductivity.grid_point_data import GridPointInput from phono3py.file_IO import read_pp_from_hdf5 from phono3py.phonon3.collision_matrix import CollisionMatrix from phono3py.phonon3.interaction import Interaction @@ -120,13 +119,13 @@ def __init__( self._pp_filename = pp_filename self._log_level = log_level - def compute(self, gp_input: GridPointInput) -> LBTECollisionResult: + def compute(self, grid_point: int) -> LBTECollisionResult: """Compute gamma and collision matrix row at one grid point. Parameters ---------- - gp_input : GridPointInput - Grid point data including BZ-grid index and band indices. + grid_point : int + BZ grid point index. Returns ------- @@ -134,7 +133,7 @@ def compute(self, gp_input: GridPointInput) -> LBTECollisionResult: gamma shape: (num_sigma, num_temp, num_band0) collision_row shape: (num_sigma, num_temp, ) """ - self._collision.set_grid_point(gp_input.grid_point) + self._collision.set_grid_point(grid_point) if self._log_level: triplets = self._pp.get_triplets_at_q()[0] @@ -143,27 +142,16 @@ def compute(self, gp_input: GridPointInput) -> LBTECollisionResult: num_sigma = len(self._sigmas) num_temp = len(self._temperatures) - num_band0 = len(gp_input.band_indices) + num_band0 = len(self._pp.band_indices) gamma_all = np.zeros((num_sigma, num_temp, num_band0), dtype="double") collision_rows: list[list[NDArray[np.double]]] = [] averaged_pp: NDArray[np.double] | None = None for i_sigma, sigma in enumerate(self._sigmas): - if self._log_level: - text = "Calculating collision matrix with " - if sigma is None: - text += "tetrahedron method." - else: - text += "sigma=%s" % sigma - if self._sigma_cutoff is None: - text += "." - else: - text += "(%4.2f SD)." % self._sigma_cutoff - print(text) self._collision.set_sigma(sigma, sigma_cutoff=self._sigma_cutoff) self._collision.run_integration_weights() - self._run_interaction(i_sigma, sigma, gp_input.grid_point) + self._run_interaction(i_sigma, sigma, grid_point) if self._is_full_pp and i_sigma == 0: averaged_pp = np.array(self._pp.averaged_interaction, dtype="double") @@ -197,9 +185,6 @@ def _run_interaction( elif i_sigma != 0 and (self._is_full_pp or self._sigma_cutoff is None): if self._log_level: print("Existing ph-ph interaction is used.") - else: - if self._log_level: - print("Calculating ph-ph interaction...") self._collision.run_interaction(is_full_pp=self._is_full_pp) def _read_interaction_from_file(self, sigma: float | None, grid_point: int) -> None: diff --git a/phono3py/conductivity/ms_smm19/velocity_providers.py b/phono3py/conductivity/ms_smm19/velocity_providers.py index f0e211fc..f4bb6c19 100644 --- a/phono3py/conductivity/ms_smm19/velocity_providers.py +++ b/phono3py/conductivity/ms_smm19/velocity_providers.py @@ -6,7 +6,7 @@ from numpy.typing import NDArray from phonopy.phonon.degeneracy import degenerate_sets -from phono3py.conductivity.grid_point_data import GridPointInput, VelocityResult +from phono3py.conductivity.grid_point_data import VelocityResult from phono3py.conductivity.ms_smm19.velocity_operator import VelocityOperator from phono3py.conductivity.utils import VOIGT_INDEX_PAIRS from phono3py.phonon.grid import ( @@ -83,14 +83,13 @@ def __init__( frequency_factor_to_THz=pp.frequency_factor_to_THz, ) - def compute(self, gp: GridPointInput) -> VelocityResult: + def compute(self, grid_point: int) -> VelocityResult: """Compute velocity operator quantities at a grid point. Parameters ---------- - gp : GridPointInput - Per-grid-point phonon data. ``frequencies`` must cover all bands - (num_band,) so that degenerate sets can be determined. + grid_point : int + BZ grid point index. Returns ------- @@ -102,16 +101,17 @@ def compute(self, gp: GridPointInput) -> VelocityResult: ``extra["velocity_operator"]`` contains the raw operator for HDF5. """ - q_point = get_qpoints_from_bz_grid_points(gp.grid_point, self._pp.bz_grid) + q_point = get_qpoints_from_bz_grid_points(grid_point, self._pp.bz_grid) self._velocity_obj.run([q_point]) assert self._velocity_obj.velocity_operators is not None # Full operator at this q: (nat3, nat3, 3) gv_op_full = self._velocity_obj.velocity_operators[0] # Filter to selected bands: (num_band0, nat3, 3) - gv_op = gv_op_full[gp.band_indices, :, :] + band_indices = self._pp.band_indices + gv_op = gv_op_full[band_indices, :, :] - gv = self._get_group_velocities(gp, gv_op_full) - vm_by_vm, kstar_order = self._get_gv_by_gv_operator(gp, gv_op) + gv = self._get_group_velocities(grid_point, gv_op_full) + vm_by_vm, kstar_order = self._get_gv_by_gv_operator(grid_point, gv_op) return VelocityResult( group_velocities=gv, vm_by_vm=vm_by_vm, @@ -124,7 +124,7 @@ def compute(self, gp: GridPointInput) -> VelocityResult: def _get_group_velocities( self, - gp: GridPointInput, + grid_point: int, gv_op_full: NDArray[np.cdouble], ) -> NDArray[np.double]: """Return diagonal group velocities with degenerate-subspace correction. @@ -138,7 +138,8 @@ def _get_group_velocities( nat3 = gv_op_full.shape[0] # Diagonal of operator gives standard group velocities gv = np.einsum("iij->ij", gv_op_full).real # (nat3, 3) - deg_sets = degenerate_sets(gp.frequencies) + frequencies = self._pp.get_phonons()[0][grid_point] + deg_sets = degenerate_sets(frequencies) pos = 0 for deg in deg_sets: if len(deg) > 1: @@ -148,18 +149,19 @@ def _get_group_velocities( gv[sl, id_dir] = np.linalg.eigvalsh(block) pos += len(deg) assert pos == nat3 - return gv[gp.band_indices, :] # (num_band0, 3) + return gv[self._pp.band_indices, :] # (num_band0, 3) def _get_gv_by_gv_operator( self, - gp: GridPointInput, + grid_point: int, gv_op: NDArray[np.cdouble], ) -> tuple[NDArray[np.cdouble], int]: """Return k-star-averaged velocity-operator outer product. Parameters ---------- - gp : GridPointInput + grid_point : int + BZ grid point index. gv_op : (num_band0, nat3, 3) complex Velocity operator filtered to selected bands. @@ -172,10 +174,10 @@ def _get_gv_by_gv_operator( """ if self._is_kappa_star: - rotation_map = get_grid_points_by_rotations(gp.grid_point, self._pp.bz_grid) + rotation_map = get_grid_points_by_rotations(grid_point, self._pp.bz_grid) else: rotation_map = get_grid_points_by_rotations( - gp.grid_point, + grid_point, self._pp.bz_grid, reciprocal_rotations=self._point_operations, ) diff --git a/phono3py/conductivity/njc23/__init__.py b/phono3py/conductivity/njc23/__init__.py index 314dbdf7..80667663 100644 --- a/phono3py/conductivity/njc23/__init__.py +++ b/phono3py/conductivity/njc23/__init__.py @@ -24,6 +24,8 @@ def _make_velocity_provider(ctx: VariantBuildContext) -> VelocityMatrixProvider: ctx.interaction, reciprocal_operations=ctx.point_operations, rotations_cartesian=ctx.rotations_cartesian, + grid_points=ctx.context.grid_points, + grid_weights=ctx.context.grid_weights, is_kappa_star=ctx.is_kappa_star, gv_delta_q=ctx.gv_delta_q, log_level=ctx.log_level, @@ -31,7 +33,7 @@ def _make_velocity_provider(ctx: VariantBuildContext) -> VelocityMatrixProvider: def _make_cv_provider(ctx: VariantBuildContext) -> HeatCapacityMatrixProvider: - return HeatCapacityMatrixProvider(ctx.interaction) + return HeatCapacityMatrixProvider(ctx.interaction, ctx.context.temperatures) def _make_rta_accumulator(ctx: VariantBuildContext) -> KuboRTAKappaAccumulator: diff --git a/phono3py/conductivity/protocols.py b/phono3py/conductivity/protocols.py index 0f82762e..0feb5052 100644 --- a/phono3py/conductivity/protocols.py +++ b/phono3py/conductivity/protocols.py @@ -18,9 +18,8 @@ Data containers --------------- -GridPointInput and provider result types (VelocityResult, -HeatCapacityResult, ScatteringResult) are defined in -``phono3py.conductivity.grid_point_data``. +Provider result types (VelocityResult, HeatCapacityResult, +ScatteringResult) are defined in ``phono3py.conductivity.grid_point_data``. """ @@ -32,7 +31,6 @@ from numpy.typing import NDArray from phono3py.conductivity.grid_point_data import ( - GridPointInput, HeatCapacityResult, ScatteringResult, VelocityResult, @@ -42,7 +40,6 @@ "VelocityProvider", "HeatCapacityProvider", "ScatteringProvider", - "GridPointInput", "VelocityResult", "HeatCapacityResult", "ScatteringResult", @@ -76,7 +73,7 @@ class VelocityProvider(Protocol): produces_gv_by_gv: bool produces_vm_by_vm: bool - def compute(self, gp: GridPointInput) -> VelocityResult: + def compute(self, grid_point: int) -> VelocityResult: """Compute velocity quantities at a grid point.""" ... @@ -94,7 +91,7 @@ class HeatCapacityProvider(Protocol): Class attributes ---------------- produces_heat_capacity_matrix : bool - True when compute_all() sets HeatCapacityResult.heat_capacity_matrix. + True when compute() sets HeatCapacityResult.heat_capacity_matrix. These flags tell the calculator which arrays to pre-allocate before the grid-point loop. @@ -103,13 +100,9 @@ class HeatCapacityProvider(Protocol): produces_heat_capacity_matrix: bool - def compute_all( + def compute( self, - frequencies: NDArray[np.double], grid_points: NDArray[np.int64], - temperatures: NDArray[np.double], - band_indices: NDArray[np.int64], - cutoff_frequency: float, ) -> HeatCapacityResult: """Compute heat-capacity quantities for all grid points at once.""" ... @@ -132,7 +125,7 @@ class ScatteringProvider(Protocol): def compute( self, - gp: GridPointInput, + grid_point: int, ) -> ScatteringResult: """Compute phonon linewidths at a grid point. diff --git a/phono3py/conductivity/rta_calculator.py b/phono3py/conductivity/rta_calculator.py index 8637c6c1..c70f3903 100644 --- a/phono3py/conductivity/rta_calculator.py +++ b/phono3py/conductivity/rta_calculator.py @@ -11,9 +11,7 @@ from phono3py.conductivity.context import ConductivityContext from phono3py.conductivity.grid_point_data import ( GridPointAggregates, - GridPointInput, compute_effective_gamma, - make_grid_point_input, ) from phono3py.conductivity.heat_capacity_providers import ModeHeatCapacityProvider from phono3py.conductivity.kappa_accumulators import RTAKappaAccumulator @@ -540,16 +538,6 @@ def _prepare_isotope_phonons(self) -> None: # Private: per-grid-point computation # ------------------------------------------------------------------ - def _make_grid_point_input(self, i_gp: int) -> GridPointInput: - return make_grid_point_input( - grid_point=int(self._context.grid_points[i_gp]), - grid_weight=int(self._context.grid_weights[i_gp]), - frequencies=self._context.frequencies, - eigenvectors=self._context.eigenvectors, - bz_grid=self._context.bz_grid, - band_indices=self._context.band_indices, - ) - def _show_log_header(self, i_gp: int) -> None: if not self._log_level: return @@ -570,13 +558,7 @@ def _show_log_header(self, i_gp: int) -> None: def _compute_bulk_heat_capacities(self) -> None: """Compute heat capacities for all grid points at once.""" assert self._context.temperatures is not None - cv_result = self._cv_provider.compute_all( - self._context.frequencies, - self._context.grid_points, - self._context.temperatures, - self._context.band_indices, - self._context.cutoff_frequency, - ) + cv_result = self._cv_provider.compute(self._context.grid_points) self._cv = cv_result.heat_capacities if cv_result.heat_capacity_matrix is not None: self._heat_capacity_matrix = cv_result.heat_capacity_matrix @@ -585,8 +567,8 @@ def _compute_all_velocities(self) -> None: """Compute velocities for all grid points.""" self._num_sampling_grid_points = 0 for i_gp in range(len(self._context.grid_points)): - gp_input = self._make_grid_point_input(i_gp) - vel_result = self._velocity_provider.compute(gp_input) + grid_point = int(self._context.grid_points[i_gp]) + vel_result = self._velocity_provider.compute(grid_point) self._num_sampling_grid_points += vel_result.num_sampling_grid_points self._gv[i_gp] = vel_result.group_velocities if vel_result.gv_by_gv is not None: @@ -598,8 +580,8 @@ def _compute_all_isotope(self) -> None: """Compute isotope scattering for all grid points.""" assert self._isotope_provider is not None for i_gp in range(len(self._context.grid_points)): - gp_input = self._make_grid_point_input(i_gp) - gamma_iso = self._isotope_provider.compute(gp_input) + grid_point = int(self._context.grid_points[i_gp]) + gamma_iso = self._isotope_provider.compute(grid_point) self._gamma_iso[:, i_gp, :] = gamma_iso[:, self._context.band_indices] def _compute_gamma_at_grid_point(self, i_gp: int) -> None: @@ -607,8 +589,8 @@ def _compute_gamma_at_grid_point(self, i_gp: int) -> None: self._show_log_header(i_gp) if not self._read_gamma: - gp_input = self._make_grid_point_input(i_gp) - scat_result = self._scattering_provider.compute(gp_input) + grid_point = int(self._context.grid_points[i_gp]) + scat_result = self._scattering_provider.compute(grid_point) gamma = scat_result.gamma ave_pp = scat_result.averaged_pp_interaction if self._is_N_U or self._is_gamma_detail: diff --git a/phono3py/conductivity/scattering_providers.py b/phono3py/conductivity/scattering_providers.py index 328afd7f..4ae13c61 100644 --- a/phono3py/conductivity/scattering_providers.py +++ b/phono3py/conductivity/scattering_providers.py @@ -8,7 +8,7 @@ from numpy.typing import NDArray from phonopy.physical_units import get_physical_units -from phono3py.conductivity.grid_point_data import GridPointInput, ScatteringResult +from phono3py.conductivity.grid_point_data import ScatteringResult from phono3py.other.isotope import Isotope from phono3py.phonon3.imag_self_energy import ImagSelfEnergy, average_by_degeneracy from phono3py.phonon3.interaction import Interaction @@ -110,13 +110,13 @@ def gamma_detail_at_q(self) -> NDArray[np.double] | None: """Return per-triplet gamma from last compute call.""" return self._gamma_detail_at_q - def compute(self, gp: GridPointInput) -> ScatteringResult: + def compute(self, grid_point: int) -> ScatteringResult: """Compute ph-ph linewidth at a grid point. Parameters ---------- - gp : GridPointInput - Per-grid-point phonon data. + grid_point : int + BZ grid point index. Returns ------- @@ -124,7 +124,7 @@ def compute(self, gp: GridPointInput) -> ScatteringResult: ``gamma`` (num_sigma, num_temp, num_band0) is set. ``averaged_pp_interaction`` (num_band0) is set when applicable. """ - num_band0 = len(gp.band_indices) + num_band0 = len(self._pp.band_indices) num_temp = len(self._temperatures) num_sigma = len(self._sigmas) @@ -139,7 +139,7 @@ def compute(self, gp: GridPointInput) -> ScatteringResult: self._averaged_pp_interaction: NDArray[np.double] | None = None - self._collision.set_grid_point(gp.grid_point) + self._collision.set_grid_point(grid_point) if self._log_level: triplets_at_q = self._pp.get_triplets_at_q()[0] @@ -147,9 +147,9 @@ def compute(self, gp: GridPointInput) -> ScatteringResult: print("Number of triplets: %d" % len(triplets_at_q), flush=True) if self._requires_full_gamma_path(): - self._run_sigmas(gp, gamma) + self._run_sigmas(grid_point, gamma) else: - self._run_sigmas_lowmem(gp, gamma) + self._run_sigmas_lowmem(grid_point, gamma) return ScatteringResult( gamma=gamma, @@ -172,24 +172,24 @@ def _requires_full_gamma_path(self) -> bool: def _run_sigmas( self, - gp: GridPointInput, + grid_point: int, gamma: NDArray[np.double], ) -> None: for j, sigma in enumerate(self._sigmas): self._collision.set_sigma(sigma, sigma_cutoff=self._sigma_cutoff) self._collision.run_integration_weights() - self._set_interaction_strength(gp, j, sigma) + self._set_interaction_strength(grid_point, j, sigma) self._allocate_gamma_detail_if_needed() - self._run_temperatures(gp, j, gamma) + self._run_temperatures(j, gamma) def _set_interaction_strength( self, - gp: GridPointInput, + grid_point: int, i_sigma: int, sigma: float | None, ) -> None: if self._read_pp: - self._set_from_file(gp, sigma) + self._set_from_file(grid_point, sigma) elif self._use_ave_pp: assert self._averaged_pp_interaction is not None self._collision.set_averaged_pp_interaction(self._averaged_pp_interaction) @@ -206,18 +206,16 @@ def _set_interaction_strength( if self._log_level: print("Existing ph-ph interaction is used.") else: - if self._log_level: - print("Calculating ph-ph interaction...") self._collision.run_interaction(is_full_pp=self._is_full_pp) if self._is_full_pp: self._averaged_pp_interaction = self._pp.averaged_interaction - def _set_from_file(self, gp: GridPointInput, sigma: float | None) -> None: + def _set_from_file(self, grid_point: int, sigma: float | None) -> None: from phono3py.file_IO import read_pp_from_hdf5 pp, _g_zero = read_pp_from_hdf5( self._pp.mesh_numbers, - grid_point=gp.grid_point, + grid_point=grid_point, sigma=sigma, sigma_cutoff=self._sigma_cutoff, filename=self._pp_filename, @@ -241,7 +239,6 @@ def _allocate_gamma_detail_if_needed(self) -> None: def _run_temperatures( self, - gp: GridPointInput, i_sigma: int, gamma: NDArray[np.double], ) -> None: @@ -263,7 +260,7 @@ def _run_temperatures( def _run_sigmas_lowmem( self, - gp: GridPointInput, + grid_point: int, gamma: NDArray[np.double], ) -> None: """Compute gamma without storing full ph-ph interaction strength.""" @@ -391,7 +388,7 @@ def _run_sigmas_lowmem( else: col = collisions - freq_at_gp = frequencies[gp.grid_point] + freq_at_gp = frequencies[grid_point] for k in range(len(self._temperatures)): gamma[j, k] = average_by_degeneracy( col[k] * col_unit_conv * pp_unit_conv, @@ -444,13 +441,13 @@ def isotope(self) -> Isotope: """Return the wrapped Isotope instance.""" return self._isotope - def compute(self, gp: GridPointInput) -> NDArray[np.double]: + def compute(self, grid_point: int) -> NDArray[np.double]: """Compute isotope linewidth at a grid point. Parameters ---------- - gp : GridPointInput - Per-grid-point phonon data. + grid_point : int + BZ grid point index. Returns ------- @@ -460,7 +457,7 @@ def compute(self, gp: GridPointInput) -> NDArray[np.double]: gamma_iso = [] for sigma in self._sigmas: self._isotope.sigma = sigma - self._isotope.set_grid_point(gp.grid_point) + self._isotope.set_grid_point(grid_point) self._isotope.run() gamma_iso.append(self._isotope.gamma) diff --git a/phono3py/conductivity/smm19/__init__.py b/phono3py/conductivity/smm19/__init__.py index f793d2e1..09ac053b 100644 --- a/phono3py/conductivity/smm19/__init__.py +++ b/phono3py/conductivity/smm19/__init__.py @@ -21,6 +21,8 @@ def _make_velocity_provider(ctx: VariantBuildContext) -> VelocityMatrixProvider: ctx.interaction, reciprocal_operations=ctx.point_operations, rotations_cartesian=ctx.rotations_cartesian, + grid_points=ctx.context.grid_points, + grid_weights=ctx.context.grid_weights, is_kappa_star=ctx.is_kappa_star, gv_delta_q=ctx.gv_delta_q, log_level=ctx.log_level, @@ -28,7 +30,7 @@ def _make_velocity_provider(ctx: VariantBuildContext) -> VelocityMatrixProvider: def _make_cv_provider(ctx: VariantBuildContext) -> ModeHeatCapacityProvider: - return ModeHeatCapacityProvider(ctx.interaction) + return ModeHeatCapacityProvider(ctx.interaction, ctx.context.temperatures) def _make_rta_accumulator(ctx: VariantBuildContext) -> SMM19RTAKappaAccumulator: diff --git a/phono3py/conductivity/velocity_providers.py b/phono3py/conductivity/velocity_providers.py index 8ee10ac8..b4fe065e 100644 --- a/phono3py/conductivity/velocity_providers.py +++ b/phono3py/conductivity/velocity_providers.py @@ -8,7 +8,7 @@ from numpy.typing import NDArray from phonopy.phonon.group_velocity import GroupVelocity -from phono3py.conductivity.grid_point_data import GridPointInput, VelocityResult +from phono3py.conductivity.grid_point_data import VelocityResult from phono3py.conductivity.utils import VOIGT_INDEX_PAIRS from phono3py.phonon.grid import ( get_grid_points_by_rotations, @@ -110,6 +110,8 @@ def __init__( pp: Interaction, point_operations: NDArray[np.int64], rotations_cartesian: NDArray[np.double], + grid_points: NDArray[np.int64], + grid_weights: NDArray[np.int64], is_kappa_star: bool = True, average_gv_over_kstar: bool = False, gv_delta_q: float | None = None, @@ -121,6 +123,7 @@ def __init__( self._pp = pp self._point_operations = point_operations self._rotations_cartesian = rotations_cartesian + self._grid_weight_map = dict(zip(grid_points, grid_weights, strict=True)) self._is_kappa_star = is_kappa_star self._average_gv_over_kstar = average_gv_over_kstar self._gv_delta_q = gv_delta_q @@ -137,13 +140,13 @@ def gv_delta_q(self) -> float | None: """Return delta-q used for finite-difference group velocity.""" return self._gv_delta_q - def compute(self, gp: GridPointInput) -> VelocityResult: + def compute(self, grid_point: int) -> VelocityResult: """Compute group velocity and v x v product at a grid point. Parameters ---------- - gp : GridPointInput - Per-grid-point phonon data. + grid_point : int + BZ grid point index. Returns ------- @@ -151,8 +154,9 @@ def compute(self, gp: GridPointInput) -> VelocityResult: ``group_velocities`` (num_band0, 3), ``gv_by_gv`` (num_band0, 6), and ``num_sampling_grid_points`` are set. """ - gv = self._get_gv(gp) - gv_by_gv, kstar_order = self._get_gv_by_gv(gp, gv) + gv = self._get_gv(grid_point) + grid_weight = int(self._grid_weight_map[grid_point]) + gv_by_gv, kstar_order = self._get_gv_by_gv(grid_point, grid_weight, gv) return VelocityResult( group_velocities=gv, gv_by_gv=gv_by_gv, @@ -163,10 +167,10 @@ def compute(self, gp: GridPointInput) -> VelocityResult: # Private helpers # ------------------------------------------------------------------ - def _get_gv(self, gp: GridPointInput) -> NDArray[np.double]: + def _get_gv(self, grid_point: int) -> NDArray[np.double]: if self._average_gv_over_kstar and len(self._point_operations) > 1: - return self._get_averaged_gv_over_kstar(gp) - return self._get_gv_at_bz_grid_point(gp.grid_point) + return self._get_averaged_gv_over_kstar(grid_point) + return self._get_gv_at_bz_grid_point(grid_point) def _get_gv_at_bz_grid_point(self, bz_gp: int) -> NDArray[np.double]: self._velocity_obj.run( @@ -175,27 +179,25 @@ def _get_gv_at_bz_grid_point(self, bz_gp: int) -> NDArray[np.double]: assert self._velocity_obj.group_velocities is not None return self._velocity_obj.group_velocities[0, self._pp.band_indices, :] - def _get_averaged_gv_over_kstar(self, gp: GridPointInput) -> NDArray[np.double]: + def _get_averaged_gv_over_kstar(self, grid_point: int) -> NDArray[np.double]: gps_rotated = get_grid_points_by_rotations( - gp.grid_point, self._pp.bz_grid, with_surface=True + grid_point, self._pp.bz_grid, with_surface=True ) assert len(gps_rotated) == len(self._point_operations) gvs = { bz_gp: self._get_gv_at_bz_grid_point(int(bz_gp)) for bz_gp in np.unique(gps_rotated) } - gv = np.zeros_like(gvs[gp.grid_point]) + gv = np.zeros_like(gvs[grid_point]) for bz_gp, r in zip(gps_rotated, self._rotations_cartesian, strict=True): gv += np.dot(gvs[bz_gp], r) return gv / len(self._point_operations) def _get_gv_by_gv( - self, gp: GridPointInput, gv: NDArray[np.double] + self, grid_point: int, grid_weight: int, gv: NDArray[np.double] ) -> tuple[NDArray[np.double], int]: if self._is_kappa_star: - multi = get_multiplicity_at_q( - gp.grid_point, self._pp, self._point_operations - ) + multi = get_multiplicity_at_q(grid_point, self._pp, self._point_operations) else: multi = 1 @@ -206,7 +208,7 @@ def _get_gv_by_gv( gv_by_gv_3x3 /= multi kstar_order = get_kstar_order( - gp.grid_weight, + grid_weight, multi, self._point_operations, verbose=self._log_level > 0, @@ -265,6 +267,8 @@ def __init__( pp: Interaction, reciprocal_operations: NDArray[np.int64], rotations_cartesian: NDArray[np.double], + grid_points: NDArray[np.int64], + grid_weights: NDArray[np.int64], is_kappa_star: bool = True, gv_delta_q: float | None = None, log_level: int = 0, @@ -274,6 +278,7 @@ def __init__( raise RuntimeError("Interaction.init_dynamical_matrix() must be called.") self._pp = pp self._is_kappa_star = is_kappa_star + self._grid_weight_map = dict(zip(grid_points, grid_weights, strict=True)) self._log_level = log_level self._velocity_obj = VelocityMatrix( pp.dynamical_matrix, @@ -285,13 +290,13 @@ def __init__( self._reciprocal_operations = reciprocal_operations self._rotations_cartesian = rotations_cartesian - def compute(self, gp: GridPointInput) -> VelocityResult: + def compute(self, grid_point: int) -> VelocityResult: """Compute velocity matrix quantities at a grid point. Parameters ---------- - gp : GridPointInput - Per-grid-point phonon data. + grid_point : int + BZ grid point index. Returns ------- @@ -302,12 +307,13 @@ def compute(self, gp: GridPointInput) -> VelocityResult: ``num_sampling_grid_points`` are set. """ - q_point = get_qpoints_from_bz_grid_points(gp.grid_point, self._pp.bz_grid) + q_point = get_qpoints_from_bz_grid_points(grid_point, self._pp.bz_grid) self._velocity_obj.run([q_point]) assert self._velocity_obj.velocity_matrices is not None gv = self._velocity_obj.group_velocities[0] + grid_weight = int(self._grid_weight_map[grid_point]) vm_by_vm, kstar_order = self._get_vm_by_vm( - gp, self._velocity_obj.velocity_matrices[0] + grid_point, grid_weight, self._velocity_obj.velocity_matrices[0] ) return VelocityResult( group_velocities=gv, @@ -321,7 +327,8 @@ def compute(self, gp: GridPointInput) -> VelocityResult: def _get_vm_by_vm( self, - gp: GridPointInput, + grid_point: int, + grid_weight: int, vm: NDArray[np.cdouble], ) -> tuple[NDArray[np.cdouble], int]: r"""Compute k-star-averaged outer product of velocity matrix elements. @@ -334,8 +341,6 @@ def _get_vm_by_vm( Returns ------- - gv : (num_band0, 3) real - Diagonal of the velocity matrix (standard group velocities). vm_by_vm : (num_band0, num_band, 6) complex k-star-averaged outer product in Voigt order (xx, yy, zz, yz, xz, xy). kstar_order : int @@ -344,7 +349,7 @@ def _get_vm_by_vm( """ if self._is_kappa_star: multi = get_multiplicity_at_q( - gp.grid_point, self._pp, self._reciprocal_operations + grid_point, self._pp, self._reciprocal_operations ) else: multi = 1 @@ -360,7 +365,7 @@ def _get_vm_by_vm( vm_by_vm /= multi kstar_order = get_kstar_order( - gp.grid_weight, + grid_weight, multi, self._reciprocal_operations, verbose=self._log_level > 0, diff --git a/phono3py/phonon/grid.py b/phono3py/phonon/grid.py index 34f386e6..57320bdb 100644 --- a/phono3py/phonon/grid.py +++ b/phono3py/phonon/grid.py @@ -995,7 +995,7 @@ def get_ir_grid_points( """ ir_grid_map = _get_ir_grid_map(bz_grid.D_diag, bz_grid.rotations, PS=bz_grid.PS) - (ir_grid_points, ir_grid_weights) = extract_ir_grid_points(ir_grid_map) + ir_grid_points, ir_grid_weights = extract_ir_grid_points(ir_grid_map) return ir_grid_points, ir_grid_weights, ir_grid_map diff --git a/test/conductivity/test_grid_point_data.py b/test/conductivity/test_grid_point_data.py index 256626c0..9de86753 100644 --- a/test/conductivity/test_grid_point_data.py +++ b/test/conductivity/test_grid_point_data.py @@ -5,7 +5,6 @@ import numpy as np from phono3py.conductivity.grid_point_data import ( - GridPointInput, HeatCapacityResult, ScatteringResult, VelocityResult, @@ -25,34 +24,6 @@ NUM_TEMP = 3 -def make_gp_input() -> GridPointInput: - """Return a minimal GridPointInput for testing.""" - return GridPointInput( - grid_point=10, - q_point=np.array([0.1, 0.2, 0.3]), - frequencies=np.ones(NUM_BAND, dtype="double"), - eigenvectors=np.eye(NUM_BAND, dtype="complex128"), - grid_weight=2, - band_indices=np.arange(NUM_BAND0, dtype="int64"), - ) - - -# --------------------------------------------------------------------------- -# GridPointInput -# --------------------------------------------------------------------------- - - -def test_grid_point_input_fields(): - """Test GridPointInput field shapes and values.""" - gp = make_gp_input() - assert gp.grid_point == 10 - assert gp.q_point.shape == (3,) - assert gp.frequencies.shape == (NUM_BAND,) - assert gp.eigenvectors.shape == (NUM_BAND, NUM_BAND) - assert gp.grid_weight == 2 - assert gp.band_indices.shape == (NUM_BAND0,) - - # --------------------------------------------------------------------------- # Protocol structural typing # --------------------------------------------------------------------------- @@ -61,27 +32,37 @@ def test_grid_point_input_fields(): class _DummyVelocityProvider: """Minimal implementation satisfying VelocityProvider protocol.""" - def compute(self, gp: GridPointInput) -> VelocityResult: + produces_gv_by_gv: bool = True + produces_vm_by_vm: bool = False + + def compute(self, grid_point: int, grid_weight: int) -> VelocityResult: return VelocityResult( group_velocities=np.zeros((NUM_BAND0, 3), dtype="double"), gv_by_gv=np.zeros((NUM_BAND0, 6), dtype="double"), - num_sampling_grid_points=gp.grid_weight, + num_sampling_grid_points=grid_weight, ) class _DummyHeatCapacityProvider: + produces_heat_capacity_matrix: bool = False + def compute( self, - gp: GridPointInput, + frequencies: np.ndarray, + grid_points: np.ndarray, temperatures: np.ndarray, + band_indices: np.ndarray, + cutoff_frequency: float, ) -> HeatCapacityResult: return HeatCapacityResult( - heat_capacities=np.zeros((len(temperatures), NUM_BAND0), dtype="double"), + heat_capacities=np.zeros( + (len(temperatures), len(grid_points), NUM_BAND0), dtype="double" + ), ) class _DummyScatteringProvider: - def compute(self, gp: GridPointInput) -> ScatteringResult: + def compute(self, grid_point: int) -> ScatteringResult: return ScatteringResult( gamma=np.zeros((1, NUM_TEMP, NUM_BAND0), dtype="double"), ) @@ -90,26 +71,26 @@ def compute(self, gp: GridPointInput) -> ScatteringResult: def test_velocity_provider_protocol_satisfied(): """A class with the right signature is accepted as VelocityProvider.""" provider: VelocityProvider = _DummyVelocityProvider() - gp = make_gp_input() - result = provider.compute(gp) + result = provider.compute(10, 2) assert result.group_velocities is not None - assert result.num_sampling_grid_points == gp.grid_weight + assert result.num_sampling_grid_points == 2 def test_heat_capacity_provider_protocol_satisfied(): """Test HeatCapacityProvider protocol with a dummy implementation.""" provider: HeatCapacityProvider = _DummyHeatCapacityProvider() - gp = make_gp_input() + freqs = np.ones((100, NUM_BAND), dtype="double") + gps = np.array([0, 1, 2], dtype="int64") temps = np.array([100.0, 200.0, 300.0]) - result = provider.compute(gp, temps) + bands = np.arange(NUM_BAND0, dtype="int64") + result = provider.compute(freqs, gps, temps, bands, 0.0) assert result.heat_capacities is not None - assert result.heat_capacities.shape == (3, NUM_BAND0) + assert result.heat_capacities.shape == (3, 3, NUM_BAND0) def test_scattering_provider_protocol_satisfied(): """Test ScatteringProvider protocol with a dummy implementation.""" provider: ScatteringProvider = _DummyScatteringProvider() - gp = make_gp_input() - result = provider.compute(gp) + result = provider.compute(10) assert result.gamma is not None assert result.gamma.shape == (1, NUM_TEMP, NUM_BAND0)