Merge pull request #208 from phonopy/kappadosthm

Add KappaDOSTHM to replace KappaDOS

phono3py/other/kaccum.py

@@ -1,6 +1,10 @@
 """Calculated accumulated property with respect to other property."""
 
-from typing import Optional
+from __future__ import annotations
+
+import warnings
+from collections.abc import Sequence
+from typing import Optional, Union
 
 import numpy as np
 from phonopy.phonon.dos import NormalDistribution
@@ -11,34 +15,57 @@
 epsilon = 1.0e-8
 
 
-class KappaDOS:
-    """Class to calculate thermal conductivity spectram."""
+class KappaDOSTHM:
+    """Class to calculate DOS like spectram with tetrahedron method.
+
+    To compute usual DOS,
+
+    kappados = KappaDOSTHM(
+        np.ones(freqs.shape)[None, :, :, None],
+        freqs,
+        bzgrid,
+        ir_grid_points,
+        ir_grid_weights=ir_weights,
+        ir_grid_map=ir_grid_map,
+        num_sampling_points=201
+    )
+
+    """
 
     def __init__(
         self,
-        mode_kappa,
-        frequencies,
+        mode_kappa: np.ndarray,
+        frequencies: np.ndarray,
         bz_grid: BZGrid,
-        ir_grid_points,
-        ir_grid_map=None,
-        frequency_points=None,
-        num_sampling_points=100,
+        ir_grid_points: np.ndarray,
+        ir_grid_weights: Optional[np.ndarray] = None,
+        ir_grid_map: Optional[np.ndarray] = None,
+        frequency_points: Optional[Union[np.ndarray, Sequence]] = None,
+        num_sampling_points: int = 100,
     ):
         """Init method.
 
+        Parameters
+        ----------
         mode_kappa : ndarray
             Target value.
             shape=(temperatures, ir_grid_points, num_band, num_elem),
             dtype='double'
         frequencies : ndarray
+            Frequencies at ir-grid points.
             shape=(ir_grid_points, num_band), dtype='double'
         bz_grid : BZGrid
         ir_grid_points : ndarray
-            Ir-grid point indices in BZ-grid.
-            shape=(ir_grid_points, ), dtype='int_'
-        ir_grid_map : ndarray, optional, default=None
-            Mapping table to ir-grid point indices in GR-grid.
-            None gives `np.arange(len(frequencies), 'int_')`.
+            Irreducible grid point indices in GR-grid.
+            shape=(num_ir_grid_points, ), dtype='int_'
+        ir_grid_weights : ndarray
+            Weights of irreducible grid points. Its sum is the number of
+            grid points in GR-grid (prod(D_diag)).
+            shape=(num_ir_grid_points, ), dtype='int_'
+        ir_grid_map : ndarray
+            Index mapping table to irreducible grid points from all grid points
+            such as, [0, 0, 2, 3, 3, ...].
+            shape=(prod(D_diag), ), dtype='int_'
         frequency_points : array_like, optional, default=None
             This is used as the frequency points. When None,
             frequency points are created from `num_sampling_points`.
@@ -64,18 +91,22 @@ def __init__(
             bzgp2irgp_map = None
         else:
             bzgp2irgp_map = self._get_bzgp2irgp_map(bz_grid, ir_grid_map)
+        if ir_grid_weights is None:
+            grid_weights = np.ones(mode_kappa.shape[1])
+        else:
+            grid_weights = ir_grid_weights
         for j, function in enumerate(("J", "I")):
             iweights = get_integration_weights(
                 self._frequency_points,
                 frequencies,
                 bz_grid,
-                grid_points=ir_grid_points,
+                grid_points=bz_grid.grg2bzg[ir_grid_points],
                 bzgp2irgp_map=bzgp2irgp_map,
                 function=function,
             )
             for i, iw in enumerate(iweights):
                 self._kdos[:, :, j] += np.transpose(
-                    np.dot(iw, mode_kappa[:, i]), axes=(1, 0, 2)
+                    np.dot(iw, mode_kappa[:, i] * grid_weights[i]), axes=(1, 0, 2)
                 )
         self._kdos /= np.prod(bz_grid.D_diag)
 
@@ -103,6 +134,37 @@ def _get_bzgp2irgp_map(self, bz_grid, ir_grid_map):
         return bzgp2irgp_map
 
 
+class KappaDOS(KappaDOSTHM):
+    """Deprecated class to calculate DOS like spectram with tetrahedron method."""
+
+    def __init__(
+        self,
+        mode_kappa: np.ndarray,
+        frequencies: np.ndarray,
+        bz_grid: BZGrid,
+        ir_grid_points: np.ndarray,
+        ir_grid_map: Optional[np.ndarray] = None,
+        frequency_points: Optional[np.ndarray] = None,
+        num_sampling_points: int = 100,
+    ):
+        """Init method."""
+        warnings.warn(
+            "KappaDOS is deprecated."
+            "Use KappaDOSTHM instead with replacing ir_grid_points in BZ-grid "
+            "by ir_grid_points in GR-grid.",
+            DeprecationWarning,
+        )
+        super().__init__(
+            mode_kappa,
+            frequencies,
+            bz_grid,
+            bz_grid.bzg2grg[ir_grid_points],
+            ir_grid_map=ir_grid_map,
+            frequency_points=frequency_points,
+            num_sampling_points=num_sampling_points,
+        )
+
+
 class GammaDOSsmearing:
     """Class to calculate Gamma spectram by smearing method."""
 

test/other/test_kaccum.py

@@ -1,13 +1,19 @@
 """Test for kaccum.py."""
 
+from __future__ import annotations
+
+from typing import Optional
+
 import numpy as np
+import pytest
 
 from phono3py import Phono3py
-from phono3py.other.kaccum import GammaDOSsmearing, KappaDOS, get_mfp
+from phono3py.other.kaccum import GammaDOSsmearing, KappaDOS, KappaDOSTHM, get_mfp
 from phono3py.phonon.grid import get_ir_grid_points
 
 
-def test_kappados_si(si_pbesol: Phono3py):
+@pytest.mark.parametrize("use_legacy_class", [False, True])
+def test_kappados_si(si_pbesol: Phono3py, use_legacy_class: bool):
     """Test KappaDOS class with Si.
 
     * 3x3 tensor vs frequency
@@ -101,7 +107,10 @@ def test_kappados_si(si_pbesol: Phono3py):
     tc = ph3.thermal_conductivity
     freq_points_in = np.array(kappados_si).reshape(-1, 3)[:, 0]
     freq_points, kdos = _calculate_kappados(
-        ph3, tc.mode_kappa[0], freq_points=freq_points_in
+        ph3,
+        tc.mode_kappa[0],
+        freq_points=freq_points_in,
+        use_legacy_class=use_legacy_class,
     )
     # for f, (jval, ival) in zip(freq_points, kdos):
     #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
@@ -110,7 +119,10 @@ def test_kappados_si(si_pbesol: Phono3py):
     )
 
     freq_points, kdos = _calculate_kappados(
-        ph3, tc.gamma[0, :, :, :, None], freq_points=freq_points_in
+        ph3,
+        tc.gamma[0, :, :, :, None],
+        freq_points=freq_points_in,
+        use_legacy_class=use_legacy_class,
     )
     # for f, (jval, ival) in zip(freq_points, kdos):
     #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
@@ -119,15 +131,18 @@ def test_kappados_si(si_pbesol: Phono3py):
     )
 
     mfp_points_in = np.array(mfpdos_si).reshape(-1, 3)[:, 0]
-    mfp_points, mfpdos = _calculate_mfpdos(ph3, mfp_points_in)
+    mfp_points, mfpdos = _calculate_mfpdos(
+        ph3, mfp_points_in, use_legacy_class=use_legacy_class
+    )
     # for f, (jval, ival) in zip(mfp_points, mfpdos):
     #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
     np.testing.assert_allclose(
         mfpdos_si, np.vstack((mfp_points, mfpdos.T)).T, rtol=0, atol=0.5
     )
 
 
-def test_kappados_nacl(nacl_pbe: Phono3py):
+@pytest.mark.parametrize("use_legacy_class", [False, True])
+def test_kappados_nacl(nacl_pbe: Phono3py, use_legacy_class: bool):
     """Test KappaDOS class with NaCl.
 
     * 3x3 tensor vs frequency
@@ -230,7 +245,10 @@ def test_kappados_nacl(nacl_pbe: Phono3py):
     )
 
     freq_points, kdos = _calculate_kappados(
-        ph3, tc.gamma[0, :, :, :, None], freq_points=freq_points_in
+        ph3,
+        tc.gamma[0, :, :, :, None],
+        freq_points=freq_points_in,
+        use_legacy_class=use_legacy_class,
     )
     for f, (jval, ival) in zip(freq_points, kdos):
         print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
@@ -239,7 +257,9 @@ def test_kappados_nacl(nacl_pbe: Phono3py):
     )
 
     mfp_points_in = np.array(mfpdos_nacl).reshape(-1, 3)[:, 0]
-    mfp_points, mfpdos = _calculate_mfpdos(ph3, mfp_points_in)
+    mfp_points, mfpdos = _calculate_mfpdos(
+        ph3, mfp_points_in, use_legacy_class=use_legacy_class
+    )
     # for f, (jval, ival) in zip(mfp_points, mfpdos):
     #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
     np.testing.assert_allclose(
@@ -281,24 +301,45 @@ def test_GammaDOSsmearing(nacl_pbe: Phono3py):
     )
 
 
-def _calculate_kappados(ph3: Phono3py, mode_prop, freq_points=None):
+def _calculate_kappados(
+    ph3: Phono3py,
+    mode_prop: np.ndarray,
+    freq_points: Optional[np.ndarray] = None,
+    use_legacy_class: bool = False,
+):
     tc = ph3.thermal_conductivity
     bz_grid = ph3.grid
     frequencies, _, _ = ph3.get_phonon_data()
-    kappados = KappaDOS(
-        mode_prop, frequencies, bz_grid, tc.grid_points, frequency_points=freq_points
-    )
+    with pytest.deprecated_call():
+        kappados = KappaDOS(
+            mode_prop,
+            frequencies,
+            bz_grid,
+            tc.grid_points,
+            frequency_points=freq_points,
+        )
     freq_points, kdos = kappados.get_kdos()
 
     ir_grid_points, _, ir_grid_map = get_ir_grid_points(bz_grid)
-    kappados = KappaDOS(
-        mode_prop,
-        tc.frequencies,
-        bz_grid,
-        tc.grid_points,
-        ir_grid_map=ir_grid_map,
-        frequency_points=freq_points,
-    )
+    if use_legacy_class:
+        with pytest.deprecated_call():
+            kappados = KappaDOS(
+                mode_prop,
+                tc.frequencies,
+                bz_grid,
+                tc.grid_points,
+                ir_grid_map=ir_grid_map,
+                frequency_points=freq_points,
+            )
+    else:
+        kappados = KappaDOSTHM(
+            mode_prop,
+            tc.frequencies,
+            bz_grid,
+            bz_grid.bzg2grg[tc.grid_points],
+            ir_grid_map=ir_grid_map,
+            frequency_points=freq_points,
+        )
     ir_freq_points, ir_kdos = kappados.get_kdos()
     np.testing.assert_equal(bz_grid.bzg2grg[tc.grid_points], ir_grid_points)
     np.testing.assert_allclose(ir_freq_points, freq_points, rtol=0, atol=1e-5)
@@ -307,20 +348,36 @@ def _calculate_kappados(ph3: Phono3py, mode_prop, freq_points=None):
     return freq_points, kdos[0, :, :, 0]
 
 
-def _calculate_mfpdos(ph3: Phono3py, mfp_points=None):
+def _calculate_mfpdos(
+    ph3: Phono3py,
+    mfp_points=None,
+    use_legacy_class: bool = False,
+):
     tc = ph3.thermal_conductivity
     bz_grid = ph3.grid
     mean_freepath = get_mfp(tc.gamma[0], tc.group_velocities)
     _, _, ir_grid_map = get_ir_grid_points(bz_grid)
-    mfpdos = KappaDOS(
-        tc.mode_kappa[0],
-        mean_freepath[0],
-        bz_grid,
-        tc.grid_points,
-        ir_grid_map=ir_grid_map,
-        frequency_points=mfp_points,
-        num_sampling_points=10,
-    )
+    if use_legacy_class:
+        with pytest.deprecated_call():
+            mfpdos = KappaDOS(
+                tc.mode_kappa[0],
+                mean_freepath[0],
+                bz_grid,
+                tc.grid_points,
+                ir_grid_map=ir_grid_map,
+                frequency_points=mfp_points,
+                num_sampling_points=10,
+            )
+    else:
+        mfpdos = KappaDOSTHM(
+            tc.mode_kappa[0],
+            mean_freepath[0],
+            bz_grid,
+            bz_grid.bzg2grg[tc.grid_points],
+            ir_grid_map=ir_grid_map,
+            frequency_points=mfp_points,
+            num_sampling_points=10,
+        )
     freq_points, kdos = mfpdos.get_kdos()
 
     return freq_points, kdos[0, :, :, 0]