Refactor kaccum

phono3py/cui/kaccum_script.py

@@ -85,7 +85,6 @@ def _get_ir_grid_info(bz_grid: BZGrid, weights, qpoints=None, ir_grid_points=Non
 
     """
     (ir_grid_points_ref, weights_ref, ir_grid_map) = get_ir_grid_points(bz_grid)
-    ir_grid_points_ref = np.array(bz_grid.grg2bzg[ir_grid_points_ref], dtype="int_")
 
     _assert_grid_in_hdf5(
         weights, qpoints, ir_grid_points, weights_ref, ir_grid_points_ref, bz_grid
@@ -300,13 +299,13 @@ def main():
     if "weight" in f_kappa:
         frequencies = f_kappa["frequency"][:]
         ir_weights = f_kappa["weight"][:]
-        ir_grid_points = None
+        ir_grid_points_BZ = None
         qpoints = f_kappa["qpoint"][:]
     elif "ir_grid_weights" in f_kappa and not args.no_gridsym:
         ir_weights = f_kappa["ir_grid_weights"][:]
-        ir_grid_points = f_kappa["ir_grid_points"][:]
+        ir_grid_points_BZ = f_kappa["ir_grid_points"][:]
         qpoints = None
-        frequencies = np.array(f_kappa["frequency"][ir_grid_points], dtype="double")
+        frequencies = np.array(f_kappa["frequency"][ir_grid_points_BZ], dtype="double")
     else:
         frequencies = f_kappa["frequency"][:]
         ir_weights = np.ones(len(frequencies), dtype="int_")
@@ -318,14 +317,14 @@ def main():
         store_dense_gp_map=True,
     )
     if args.no_gridsym or (ir_weights == 1).all():
-        ir_grid_points = bz_grid.grg2bzg
+        ir_grid_points = None
         ir_grid_map = None
     else:
         ir_grid_points, ir_grid_map = _get_ir_grid_info(
             bz_grid,
             ir_weights,
             qpoints=qpoints,
-            ir_grid_points=ir_grid_points,
+            ir_grid_points=bz_grid.bzg2grg[ir_grid_points_BZ],
         )
     conditions = frequencies > 0
     if np.logical_not(conditions).sum() > 3:

phono3py/other/kaccum.py

@@ -20,12 +20,9 @@ class KappaDOSTHM:
     To compute usual DOS,
 
     kappados = KappaDOSTHM(
-        np.ones(freqs.shape)[None, :, :, None],
+        np.ones(freqs.shape, dtype=float)[None, :, :, None],
         freqs,
         bzgrid,
-        ir_grid_points,
-        ir_grid_weights=ir_weights,
-        ir_grid_map=ir_grid_map,
         num_sampling_points=201
     )
 
@@ -36,7 +33,7 @@ class KappaDOSTHM:
         mode_kappa: np.ndarray,
         frequencies: np.ndarray,
         bz_grid: BZGrid,
-        ir_grid_points: np.ndarray,
+        ir_grid_points: Optional[np.ndarray] = None,
         ir_grid_weights: Optional[np.ndarray] = None,
         ir_grid_map: Optional[np.ndarray] = None,
         frequency_points: Optional[Union[np.ndarray, Sequence]] = None,
@@ -54,20 +51,24 @@ class KappaDOSTHM:
             Frequencies at ir-grid points.
             shape=(ir_grid_points, num_band), dtype='double'
         bz_grid : BZGrid
+            BZGrid instance.
         ir_grid_points : ndarray
-            Irreducible grid point indices in GR-grid.
+            Irreducible grid point indices in GR-grid (as obtained by
+            get_ir_grid_points).
             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)).
+            Weights of irreducible grid points. Its sum is the number of grid
+            points in GR-grid (prod(D_diag)) (as obtained by
+            get_ir_grid_points).
             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, ...].
+            in GR-grid such as, [0, 0, 2, 3, 3, ...]. (as obtained by
+            get_ir_grid_points).
             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`.
+            This is used as the frequency points. When None, frequency points
+            are created from `num_sampling_points`.
         num_sampling_points : int, optional, default=100
             Number of uniform sampling points.
 
@@ -85,13 +86,20 @@ class KappaDOSTHM:
         self._kdos = np.zeros(
             (n_temp, len(self._frequency_points), 2, n_elem), dtype="double"
         )
-
-        if ir_grid_map is None:
-            bzgp2irgp_map = None
+        if ir_grid_points is None:
+            _ir_grid_points = np.arange(len(frequencies), dtype=int)
         else:
-            bzgp2irgp_map = self._get_bzgp2irgp_map(bz_grid, ir_grid_map)
+            _ir_grid_points = ir_grid_points
+        grid_points = bz_grid.grg2bzg[_ir_grid_points]
+        if ir_grid_map is None:
+            _ir_grid_map = np.arange(len(frequencies), dtype=int)
+        else:
+            _ir_grid_map = ir_grid_map
+        bzgp2irgp_map = self._get_bzgp2irgp_map(
+            bz_grid.bzg2grg, _ir_grid_map, _ir_grid_points
+        )
         if ir_grid_weights is None:
-            grid_weights = np.ones(mode_kappa.shape[1])
+            grid_weights = np.ones(mode_kappa.shape[1], dtype="int_")
         else:
             grid_weights = ir_grid_weights
         for j, function in enumerate(("J", "I")):
@@ -99,7 +107,7 @@ class KappaDOSTHM:
                 self._frequency_points,
                 frequencies,
                 bz_grid,
-                grid_points=bz_grid.grg2bzg[ir_grid_points],
+                grid_points=grid_points,
                 bzgp2irgp_map=bzgp2irgp_map,
                 function=function,
             )
@@ -124,11 +132,38 @@ class KappaDOSTHM:
         """
         return self._frequency_points, self._kdos
 
-    def _get_bzgp2irgp_map(self, bz_grid, ir_grid_map):
+    def _get_bzgp2irgp_map(self, bzg2grg, ir_grid_map, ir_grid_points):
+        """Return mapping table from BZ-grid indices to ir-grid point indices.
+
+        More precisely, return mapping table from grid points in BZ-grid to
+        indices of ir-grid points. The length of the set of the indices is the
+        number of the ir-grid points.
+
+        Parameters
+        ----------
+        bzg2grg : ndarray
+            Mapping table from BZ-grid to GR-grid.
+            shape=(len(all-BZ-grid-points), ), dtype='int_'
+        ir_grid_map : ndarray
+            Mapping table from all grid points to ir-grid points in GR-grid.
+            shape=(np.prod(D_diag), ), dtype='int_'
+        ir_grid_points : ndarray
+            Irreducible grid points in GR-grid. shape=(num_ir_grid_points, ),
+            dtype='int_'
+
+        Returns
+        -------
+        np.ndarray
+            Mapping table from BZ-grid to indices of ir-grid points.
+            shape=(len(ir-grid-points), ), dtype='
+
+        """
         unique_gps = np.unique(ir_grid_map)
+        assert np.array_equal(unique_gps, ir_grid_points)
+        # ir-grid points in GR-grid to the index of unique grid points.
         gp_map = {j: i for i, j in enumerate(unique_gps)}
         bzgp2irgp_map = np.array(
-            [gp_map[ir_grid_map[grgp]] for grgp in bz_grid.bzg2grg], dtype="int_"
+            [gp_map[ir_grid_map[grgp]] for grgp in bzg2grg], dtype="int_"
         )
         return bzgp2irgp_map
 
@@ -236,7 +271,7 @@ def run_prop_dos(
         mode_prop,
         frequencies,
         bz_grid,
-        bz_grid.bzg2grg[ir_grid_points],
+        ir_grid_points=ir_grid_points,
         ir_grid_map=ir_grid_map,
         num_sampling_points=num_sampling_points,
     )
@@ -266,7 +301,7 @@ def run_mfp_dos(
             mode_prop[i : i + 1, :, :],
             mean_freepath[i],
             bz_grid,
-            bz_grid.bzg2grg[ir_grid_points],
+            ir_grid_points=ir_grid_points,
             ir_grid_map=ir_grid_map,
             num_sampling_points=num_sampling_points,
         )

test/other/test_kaccum.py

@@ -323,51 +323,51 @@ def test_run_prop_dos(si_pbesol: Phono3py):
     ref_kdos = [
         0.00000,
         0.00000,
-        2.19162,
-        5.16697,
-        28.22125,
-        18.97280,
-        58.56343,
-        12.19206,
-        69.05896,
-        3.47035,
-        73.17626,
-        1.48915,
-        74.74544,
-        0.43485,
-        75.87064,
-        1.74135,
-        79.08179,
-        2.30428,
-        81.21678,
+        2.19027,
+        5.16378,
+        25.73771,
+        15.50548,
+        56.63845,
+        19.48634,
+        68.64061,
+        3.16016,
+        72.58422,
+        1.62370,
+        74.53304,
+        0.78480,
+        76.94719,
+        5.37456,
+        80.63569,
+        0.61517,
+        81.15021,
         0.00000,
     ]
     # print(",".join([f"{v:10.5f}" for v in mfp[0, :, :, 0].ravel()]))
     ref_mfp = [
         0.00000,
         0.00000,
-        29.19150,
-        0.02604,
-        42.80717,
-        0.01202,
-        52.09457,
-        0.01158,
-        61.79908,
-        0.01140,
-        69.49177,
-        0.00784,
-        74.57499,
-        0.00501,
-        77.99145,
-        0.00364,
-        80.33477,
-        0.00210,
-        81.21678,
+        33.66617,
+        0.02216,
+        44.81836,
+        0.01032,
+        53.09252,
+        0.01062,
+        62.18077,
+        0.01085,
+        69.59018,
+        0.00765,
+        74.58994,
+        0.00496,
+        77.95139,
+        0.00359,
+        80.27019,
+        0.00209,
+        81.15021,
         0.00000,
     ]
     # print(",".join([f"{v:10.5f}" for v in sampling_points[0]]))
     ref_sampling_points = [
-        -0.00000,
+        0.00000,
         1.69664,
         3.39328,
         5.08992,
@@ -381,15 +381,15 @@ def test_run_prop_dos(si_pbesol: Phono3py):
     # print(",".join([f"{v:10.5f}" for v in sampling_points_mfp[0]]))
     ref_sampling_points_mfp = [
         0.00000,
-        803.91710,
-        1607.83420,
-        2411.75130,
-        3215.66841,
-        4019.58551,
-        4823.50261,
-        5627.41971,
-        6431.33681,
-        7235.25391,
+        803.28312,
+        1606.56625,
+        2409.84937,
+        3213.13249,
+        4016.41562,
+        4819.69874,
+        5622.98186,
+        6426.26499,
+        7229.54811,
     ]
     np.testing.assert_allclose(ref_kdos, kdos[0, :, :, 0].ravel(), atol=0.5)
     np.testing.assert_allclose(ref_mfp, mfp[0, :, :, 0].ravel(), atol=0.5)