Merge pull request #206 from phonopy/boundary-mfp

Write boundary-mfp in kappa-*.hdf5 file

doc/changelog.md

@@ -2,6 +2,10 @@
 
 # Change Log
 
+## Feb-2-2024: Version 2.9.2
+
+- `boundary_mfp` value is stored in `kappa-*.hdf5` file when it is specified.
+
 ## Dec-26-2023: Version 2.9.1
 
 - Release to build conda-forge package.

doc/command-options.md

@@ -561,6 +561,8 @@ Then running with isotope calculation:
 In the result hdf5 file, currently isotope scattering strength is not written
 out, i.e., `gamma` is still imaginary part of self energy of ph-ph scattering.
 
+(boundary_mfp_option)=
+
 ### `--boundary-mfp`, `--bmfp` (`BOUNDARY_MFP`)
 
 A most simple phonon boundary scattering treatment is included. $v_g/L$ is just

doc/conf.py

@@ -60,7 +60,7 @@
 # The short X.Y version.
 version = "2.9"
 # The full version, including alpha/beta/rc tags.
-release = "2.9.1"
+release = "2.9.2"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

doc/hdf5_howto.md

@@ -123,8 +123,9 @@ The array shape is (irreducible q-point, phonon band).
 (kappa_hdf5_file_gamma)=
 ### gamma
 
-Imaginary part of self energy. The physical unit is THz, where THz
-is in the ordinal frequency not the angular frequency.
+Imaginary part of self energy of phonon bubble diagram (phonon-phonon
+scattering). The physical unit is THz, where THz is in the ordinal frequency not
+the angular frequency.
 
 The array shape for all grid-points (irreducible q-points) is
 (temperature, irreducible q-point, phonon band).
@@ -262,9 +263,16 @@ P_{\mathbf{q}j} = \frac{1}{(3n_\mathrm{a})^2} \sum_{\lambda'\lambda''}
 $$
 
 This is not going to be calculated in the RTA thermal coductivity
-calculation mode by default. To calculate this, `--full_pp` option
+calculation mode by default. To calculate this, `--full-pp` option
 has to be specified (see {ref}`full_pp_option`).
 
+### boundary_mfp
+
+A value specified by {ref}`boundary_mfp_option`. The physical unit is
+micrometre.
+
+When `--boundary-mfp` option is explicitly specified, its value is stored here.
+
 ### kappa_unit_conversion
 
 This is used to convert the physical unit of lattice thermal

phono3py/conductivity/base.py

@@ -376,6 +376,8 @@ def __init__(
         volume = self._pp.primitive.volume
         self._conversion_factor = unit_to_WmK / volume
 
+        self._averaged_pp_interaction = None
+
         # `self._velocity_obj` is the instance of an inherited class of
         # `GroupVelocity`. `self._init_velocity()` is the method setup the instance,
         # which must be implmented in the inherited class of `ConductivityBase`.
@@ -633,6 +635,11 @@ def get_averaged_pp_interaction(self):
         )
         return self.averaged_pp_interaction
 
+    @property
+    def boundary_mfp(self) -> float:
+        """Return boundary MFP."""
+        return self._boundary_mfp
+
     def get_number_of_sampling_grid_points(self):
         """Return number of grid points.
 

phono3py/conductivity/direct_solution.py

@@ -109,7 +109,6 @@ def __init__(
         self._init_velocity(gv_delta_q)
 
         self._lang = lang
-        self._averaged_pp_interaction = None
         self._collision_eigenvalues = None
         self._is_reducible_collision_matrix = is_reducible_collision_matrix
         self._solve_collective_phonon = solve_collective_phonon

phono3py/conductivity/rta.py

@@ -98,7 +98,6 @@ def __init__(
         self._gamma_detail_at_q = None
         self._use_ave_pp = use_ave_pp
         self._use_const_ave_pp = None
-        self._averaged_pp_interaction = None
         self._num_ignored_phonon_modes = None
 
         super().__init__(

phono3py/conductivity/utils.py

@@ -269,6 +269,7 @@ def write_kappa(
         qpoints = br.qpoints
         grid_points = br.grid_points
         weights = br.grid_weights
+        boundary_mfp = br.boundary_mfp
 
         for i, sigma in enumerate(sigmas):
             if kappa is None:
@@ -315,6 +316,7 @@ def write_kappa(
             write_kappa_to_hdf5(
                 temperatures,
                 mesh,
+                boundary_mfp=boundary_mfp,
                 bz_grid=bz_grid,
                 frequency=frequencies,
                 group_velocity=gv,
@@ -488,13 +490,13 @@ def write_collision(
     @staticmethod
     def write_kappa(
         lbte: "cond_LBTE_type",
-        volume,
-        is_reducible_collision_matrix=False,
-        write_LBTE_solution=False,
-        pinv_solver=None,
-        compression="gzip",
-        filename=None,
-        log_level=0,
+        volume: float,
+        is_reducible_collision_matrix: bool = False,
+        write_LBTE_solution: bool = False,
+        pinv_solver: Optional[int] = None,
+        compression: str = "gzip",
+        filename: Optional[str] = None,
+        log_level: int = 0,
     ):
         """Write kappa related properties into a hdf5 file."""
         from phono3py.conductivity.direct_solution import (
@@ -547,6 +549,7 @@ def write_kappa(
         f_vector = lbte.get_f_vectors()
         mode_cv = lbte.mode_heat_capacities
         mfp = lbte.get_mean_free_path()
+        boundary_mfp = lbte.boundary_mfp
 
         coleigs = lbte.get_collision_eigenvalues()
         # After kappa calculation, the variable is overwritten by unitary matrix
@@ -613,6 +616,7 @@ def write_kappa(
             write_kappa_to_hdf5(
                 temperatures,
                 mesh,
+                boundary_mfp=boundary_mfp,
                 bz_grid=bz_grid,
                 frequency=frequencies,
                 group_velocity=gv,

phono3py/file_IO.py

@@ -33,9 +33,12 @@
 # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
+from __future__ import annotations
 
 import os
 import warnings
+from collections.abc import Sequence
+from typing import Optional, Union
 
 import h5py
 import numpy as np
@@ -919,6 +922,7 @@ def write_collision_eigenvalues_to_hdf5(
 def write_kappa_to_hdf5(
     temperature,
     mesh,
+    boundary_mfp: float = None,
     bz_grid=None,
     frequency=None,
     group_velocity=None,
@@ -1059,6 +1063,8 @@ def write_kappa_to_hdf5(
             w.create_dataset("sigma_cutoff_width", data=sigma_cutoff)
         if kappa_unit_conversion is not None:
             w.create_dataset("kappa_unit_conversion", data=kappa_unit_conversion)
+        if boundary_mfp is not None:
+            w.create_dataset("boundary_mfp", data=boundary_mfp)
 
         if verbose:
             text = "Thermal conductivity related properties "
@@ -1721,12 +1727,12 @@ def get_length_of_first_line(f):
 
 def _get_filename_suffix(
     mesh,
-    grid_point=None,
-    band_indices=None,
-    sigma=None,
-    sigma_cutoff=None,
-    temperature=None,
-    filename=None,
+    grid_point: Optional[int] = None,
+    band_indices: Optional[Union[np.ndarray, Sequence]] = None,
+    sigma: Optional[float] = None,
+    sigma_cutoff: Optional[float] = None,
+    temperature: Optional[float] = None,
+    filename: Optional[str] = None,
 ):
     """Return filename suffix corresponding to parameters."""
     suffix = "-m%d%d%d" % tuple(mesh)

phono3py/version.py

@@ -34,4 +34,4 @@
 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
 
-__version__ = "2.9.1"
+__version__ = "2.9.2"

test/file_IO/test_file_IO.py

@@ -0,0 +1,95 @@
+"""Tests of Phono3py API."""
+
+import pathlib
+from collections.abc import Sequence
+from typing import Optional
+
+import h5py
+import numpy as np
+
+from phono3py import Phono3py
+from phono3py.file_IO import _get_filename_suffix
+
+cwd = pathlib.Path(__file__).parent
+cwd_called = pathlib.Path.cwd()
+
+
+def test_kappa_filename():
+    """Test _get_filename_suffix."""
+    mesh = [4, 4, 4]
+    grid_point = None
+    band_indices = None
+    sigma = None
+    sigma_cutoff = None
+    filename = None
+    suffix = _get_filename_suffix(
+        mesh,
+        grid_point=grid_point,
+        band_indices=band_indices,
+        sigma=sigma,
+        sigma_cutoff=sigma_cutoff,
+        filename=filename,
+    )
+    full_filename = "kappa" + suffix + ".hdf5"
+    assert full_filename == "kappa-m444.hdf5"
+
+
+def test_kappa_hdf5_with_boundary_mpf(si_pbesol: Phono3py):
+    """Test boundary_mfp in kappa-*.hdf5.
+
+    Remember to clean files created by
+    Phono3py.run_thermal_conductivity(write_kappa=True).
+
+    """
+    key_ref = [
+        "boundary_mfp",
+        "frequency",
+        "gamma",
+        "grid_point",
+        "group_velocity",
+        "gv_by_gv",
+        "heat_capacity",
+        "kappa",
+        "kappa_unit_conversion",
+        "mesh",
+        "mode_kappa",
+        "qpoint",
+        "temperature",
+        "version",
+        "weight",
+    ]
+
+    boundary_mfp = 10000.0
+    kappa_filename = _set_kappa(
+        si_pbesol, [4, 4, 4], write_kappa=True, boundary_mfp=boundary_mfp
+    )
+    file_path = pathlib.Path(cwd_called / kappa_filename)
+    with h5py.File(file_path, "r") as f:
+        np.testing.assert_almost_equal(f["boundary_mfp"][()], boundary_mfp)
+        assert set(list(f)) == set(key_ref)
+
+    if file_path.exists():
+        file_path.unlink()
+
+
+def _set_kappa(
+    ph3: Phono3py,
+    mesh: Sequence,
+    is_isotope: bool = False,
+    is_full_pp: bool = False,
+    write_kappa: bool = False,
+    boundary_mfp: Optional[float] = None,
+) -> str:
+    ph3.mesh_numbers = mesh
+    ph3.init_phph_interaction()
+    ph3.run_thermal_conductivity(
+        temperatures=[
+            300,
+        ],
+        is_isotope=is_isotope,
+        is_full_pp=is_full_pp,
+        write_kappa=write_kappa,
+        boundary_mfp=boundary_mfp,
+    )
+    suffix = _get_filename_suffix(mesh)
+    return "kappa" + suffix + ".hdf5"