Refactor kaccum

doc/auxiliary-tools.md

@@ -80,36 +80,17 @@ With $19\times 19\times 19$ mesh:
 :width: 25%
 ```
 
-###General options
-
-#### `--pa`
-
-See {ref}`pa_option`.
-
-#### `-c`
-
-Unit cell filename is specified with this option, e.g., `-c POSCAR-unitcell`.
-
-#### `--qe`
-
-Let `phono3py-kaccum` read a QE (pw) unit cell file with `-c` option, for
-example:
-
-```bash
-% phono3py-kaccum --qe kappa-m191919.hdf5
-```
+That calculated by QE with $19\times 19\times 19$ mesh:
 
 ```{image} Si-kaccum-pwscf.png
 :width: 25%
 ```
 
-#### `--crystal`
+###General options
 
-Analogous to `--qe`, but to be used with the CRYSTAL interface.
+#### `--pa`
 
-#### `--turbomole`
-
-Analogous to `--qe`, but to be used with the TURBOMOLE interface
+See {ref}`pa_option`.
 
 #### `--temperature`
 

doc/changelog.md

@@ -6,9 +6,9 @@
 
 This is a major version release. There are backward-incompatible changes.
 
-- Method to transform supercell third-order force constants fc3 in real to
-  reciprocal space was changed as described at Version 2.9.0 changelog below.
-  This results in the change of results with respect to those obtained by
+- Calculation method to transform supercell third-order force constants fc3 in
+  real to reciprocal space was changed as described at Version 2.9.0 changelog
+  below. This results in the change of results with respect to those obtained by
   phono3py version 2. To emulate v2 behaviour, use `--v2` option in phono3py
   command line script. For `Phono3py` class , `make_r0_average=True` (default)
   when instantiating it, and similarly for `phono3py.load` function.

phono3py/cui/kaccum_script.py

@@ -129,14 +129,14 @@ def _assert_grid_in_hdf5(
 def _get_calculator(args):
     """Return calculator name."""
     interface_mode = None
-    if args.qe_mode:
-        interface_mode = "qe"
-    elif args.crystal_mode:
-        interface_mode = "crystal"
-    elif args.abinit_mode:
-        interface_mode = "abinit"
-    elif args.turbomole_mode:
-        interface_mode = "turbomole"
+    # if args.qe_mode:
+    #     interface_mode = "qe"
+    # elif args.crystal_mode:
+    #     interface_mode = "crystal"
+    # elif args.abinit_mode:
+    #     interface_mode = "abinit"
+    # elif args.turbomole_mode:
+    #     interface_mode = "turbomole"
     return interface_mode
 
 
@@ -157,16 +157,6 @@ def _read_files(args):
     return cell, f
 
 
-def _read_files_by_collect_cell_info(cell_filename, interface_mode):
-    cell_info = collect_cell_info(
-        interface_mode=interface_mode,
-        cell_filename=cell_filename,
-        supercell_matrix=np.eye(3, dtype=int),
-        phonopy_yaml_cls=Phono3pyYaml,
-    )
-    return cell_info
-
-
 def _get_mode_property(args, f_kappa):
     """Read property data from hdf5 file object."""
     if args.pqj:
@@ -206,14 +196,14 @@ def _get_parser():
         default=None,
         help="Same as PRIMITIVE_AXES tags",
     )
-    parser.add_argument(
-        "-c",
-        "--cell",
-        dest="cell_filename",
-        metavar="FILE",
-        default=None,
-        help="Read unit cell",
-    )
+    # parser.add_argument(
+    #     "-c",
+    #     "--cell",
+    #     dest="cell_filename",
+    #     metavar="FILE",
+    #     default=None,
+    #     help="Read unit cell",
+    # )
     parser.add_argument(
         "--gv", action="store_true", help="Calculate for gv_x_gv (tensor)"
     )
@@ -271,24 +261,6 @@ def _get_parser():
         action="store_true",
         help="Use smearing method (only for scalar density)",
     )
-    parser.add_argument(
-        "--qe", "--pwscf", dest="qe_mode", action="store_true", help="Invoke Pwscf mode"
-    )
-    parser.add_argument(
-        "--crystal",
-        dest="crystal_mode",
-        action="store_true",
-        help="Invoke CRYSTAL mode",
-    )
-    parser.add_argument(
-        "--abinit", dest="abinit_mode", action="store_true", help="Invoke Abinit mode"
-    )
-    parser.add_argument(
-        "--turbomole",
-        dest="turbomole_mode",
-        action="store_true",
-        help="Invoke TURBOMOLE mode",
-    )
     parser.add_argument(
         "--no-gridsym",
         dest="no_gridsym",
@@ -363,8 +335,12 @@ def main():
             'Use of "phono3py-kaccum CRYSTAL_STRUCTURE_FILE" is not supported.'
         )
     else:
-        interface_mode = _get_calculator(args)
-        cell_info = _read_files_by_collect_cell_info(args.cell_filename, interface_mode)
+        cell_info = collect_cell_info(
+            supercell_matrix=np.eye(3, dtype=int),
+            phonopy_yaml_cls=Phono3pyYaml,
+        )
+        cell_filename = cell_info["optional_structure_info"][0]
+        print(f'# Crystal structure was read from "{cell_filename}".')
         cell = cell_info["unitcell"]
         phpy_yaml = cell_info.get("phonopy_yaml", None)
         if phpy_yaml is not None:

phono3py/cui/phono3py_argparse.py

@@ -282,7 +282,7 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
             help="Read third order force constants",
         )
     parser.add_argument(
-        "--v3",
+        "--v2",
         dest="is_fc3_r0_average",
         action="store_true",
         default=False,

phono3py/other/kaccum.py

@@ -209,15 +209,34 @@ def run_prop_dos(
     mode_prop,
     ir_grid_map,
     ir_grid_points,
-    num_sampling_points,
+    num_sampling_points: int,
     bz_grid: BZGrid,
 ):
-    """Run DOS-like calculation."""
+    """Run DOS-like calculation.
+
+    This is a simple wrapper of KappsDOSTHM.
+
+    Parameters
+    ----------
+    frequencies:
+        Frequencies at ir-grid points.
+    mode_prop:
+        Properties at  ir-grid points.
+    ir_grid_map:
+        Obtained by get_ir_grid_points(bz_grid)[2].
+    ir_grid_points:
+        Obtained by get_ir_grid_points(bz_grid)[0].
+    num_sampling_points:
+        Number of sampling points in horizontal axis.
+    bz_grid:
+        BZ grid.
+
+    """
     kappa_dos = KappaDOSTHM(
         mode_prop,
         frequencies,
         bz_grid,
-        ir_grid_points,
+        bz_grid.bzg2grg[ir_grid_points],
         ir_grid_map=ir_grid_map,
         num_sampling_points=num_sampling_points,
     )
@@ -227,7 +246,12 @@ def run_prop_dos(
 
 
 def run_mfp_dos(
-    mean_freepath, mode_prop, ir_grid_map, ir_grid_points, num_sampling_points, bz_grid
+    mean_freepath,
+    mode_prop,
+    ir_grid_map,
+    ir_grid_points,
+    num_sampling_points: int,
+    bz_grid: BZGrid,
 ):
     """Run DOS-like calculation for mean free path.
 
@@ -242,7 +266,7 @@ def run_mfp_dos(
             mode_prop[i : i + 1, :, :],
             mean_freepath[i],
             bz_grid,
-            ir_grid_points,
+            bz_grid.bzg2grg[ir_grid_points],
             ir_grid_map=ir_grid_map,
             num_sampling_points=num_sampling_points,
         )

test/other/test_kaccum.py

@@ -7,7 +7,13 @@ from typing import Optional
 import numpy as np
 
 from phono3py import Phono3py
-from phono3py.other.kaccum import GammaDOSsmearing, KappaDOSTHM, get_mfp
+from phono3py.other.kaccum import (
+    GammaDOSsmearing,
+    KappaDOSTHM,
+    get_mfp,
+    run_mfp_dos,
+    run_prop_dos,
+)
 from phono3py.phonon.grid import get_ir_grid_points
 
 
@@ -291,6 +297,105 @@ def test_GammaDOSsmearing(nacl_pbe: Phono3py):
     )
 
 
+def test_run_prop_dos(si_pbesol: Phono3py):
+    ph3 = si_pbesol
+    ph3.mesh_numbers = [7, 7, 7]
+    ph3.init_phph_interaction()
+    ph3.run_thermal_conductivity(
+        temperatures=[
+            300,
+        ]
+    )
+    bz_grid = ph3.grid
+    ir_grid_points, _, ir_grid_map = get_ir_grid_points(bz_grid)
+    tc = ph3.thermal_conductivity
+
+    kdos, sampling_points = run_prop_dos(
+        tc.frequencies, tc.mode_kappa[0], ir_grid_map, ir_grid_points, 10, bz_grid
+    )
+    mean_freepath = get_mfp(tc.gamma[0], tc.group_velocities)
+    mfp, sampling_points_mfp = run_mfp_dos(
+        mean_freepath, tc.mode_kappa[0], ir_grid_map, ir_grid_points, 10, bz_grid
+    )
+
+    # print(",".join([f"{v:10.5f}" for v in kdos[0, :, :, 0].ravel()]))
+    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,
+        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,
+        0.00000,
+    ]
+    # print(",".join([f"{v:10.5f}" for v in sampling_points[0]]))
+    ref_sampling_points = [
+        -0.00000,
+        1.69664,
+        3.39328,
+        5.08992,
+        6.78656,
+        8.48320,
+        10.17984,
+        11.87648,
+        13.57312,
+        15.26976,
+    ]
+    # 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,
+    ]
+    np.testing.assert_allclose(ref_kdos, kdos[0, :, :, 0].ravel(), atol=1e-2)
+    np.testing.assert_allclose(ref_mfp, mfp[0, :, :, 0].ravel(), atol=1e-2)
+    np.testing.assert_allclose(ref_sampling_points, sampling_points[0], atol=1e-4)
+    np.testing.assert_allclose(ref_sampling_points_mfp, sampling_points_mfp[0], rtol=10)
+
+
 def _calculate_kappados(
     ph3: Phono3py,
     mode_prop: np.ndarray,