From d0b7ce08be039543e6a90cd8c8d0991adf8eadab Mon Sep 17 00:00:00 2001
From: Atsushi Togo <atz.togo@gmail.com>
Date: Mon, 9 Dec 2024 15:29:52 +0900
Subject: [PATCH 1/4] Add velph-phonopy command

---
 doc/velph.md                                  |   2 +-
 src/phelel/velph/cli/__init__.py              |   5 +-
 src/phelel/velph/cli/init/__init__.py         |  16 ++
 src/phelel/velph/cli/init/init.py             | 242 ++++++++++++------
 src/phelel/velph/cli/phono3py/__init__.py     |  25 +-
 src/phelel/velph/cli/phono3py/init.py         |  22 +-
 src/phelel/velph/cli/utils.py                 |   1 +
 .../cli/phono3py/init/test_phono3py_init.py   |  30 +--
 8 files changed, 206 insertions(+), 137 deletions(-)

diff --git a/doc/velph.md b/doc/velph.md
index a36ec69..ca8989e 100644
--- a/doc/velph.md
+++ b/doc/velph.md
@@ -208,7 +208,7 @@ source /opt/intel/oneapi/setvars.sh --config="/home/togo/.oneapi-config"
 pe = "paris 24"
 
 ...
-[vasp.phono3py.phonon.scheduler]
+[vasp.phonopy.scheduler]
 scheduler_template = '''#!/bin/bash
 #QSUB2 core 192
 #QSUB2 mpi 192
diff --git a/src/phelel/velph/cli/__init__.py b/src/phelel/velph/cli/__init__.py
index 37b759c..34292f8 100644
--- a/src/phelel/velph/cli/__init__.py
+++ b/src/phelel/velph/cli/__init__.py
@@ -14,10 +14,11 @@ def cmd_root():
 from phelel.velph.cli.generate import cmd_generate  # noqa F401
 from phelel.velph.cli.hints import cmd_hints  # noqa F401
 from phelel.velph.cli.init import cmd_init  # noqa F401
-from phelel.velph.cli.phono3py import cmd_phono3py  # noqa F401
 from phelel.velph.cli.nac import cmd_nac  # noqa F401
+from phelel.velph.cli.phelel import cmd_phelel  # noqa F401
+from phelel.velph.cli.phonopy import cmd_phonopy  # noqa F401
+from phelel.velph.cli.phono3py import cmd_phono3py  # noqa F401
 from phelel.velph.cli.ph_bands import cmd_ph_bands  # noqa F401
 from phelel.velph.cli.relax import cmd_relax  # noqa F401
 from phelel.velph.cli.selfenergy import cmd_selfenergy  # noqa F401
-from phelel.velph.cli.phelel import cmd_phelel  # noqa F401
 from phelel.velph.cli.transport import cmd_transport  # noqa F401
diff --git a/src/phelel/velph/cli/init/__init__.py b/src/phelel/velph/cli/init/__init__.py
index 45c1a9f..8c1d4e9 100644
--- a/src/phelel/velph/cli/init/__init__.py
+++ b/src/phelel/velph/cli/init/__init__.py
@@ -145,6 +145,20 @@
         f"(phelel_nosym: bool, default={VelphInitParams.phelel_nosym})"
     ),
 )
+@click.option(
+    "--phonopy-max-num-atoms",
+    "phonopy_max_num_atoms",
+    nargs=1,
+    default=None,
+    type=int,
+    help=(
+        "Determine phonopy supercell dimension so that number of atoms in supercell "
+        "for phonopy is less than this number if different dimension from "
+        "that of electron-phonon (phelel) is expected. "
+        "(phonopy_max_num_atoms: int, "
+        f"default={VelphInitParams.phono3py_max_num_atoms})"
+    ),
+)
 @click.option(
     "--phono3py-max-num-atoms",
     "phono3py_max_num_atoms",
@@ -246,6 +260,7 @@ def cmd_init(
     max_num_atoms: Optional[int],
     phelel_dir_name: str,
     phelel_nosym: Optional[bool],
+    phonopy_max_num_atoms: Optional[int],
     phono3py_max_num_atoms: Optional[int],
     primitive_cell_choice: Optional[str],
     project_folder: str,
@@ -287,6 +302,7 @@ def cmd_init(
         "magmom": magmom,
         "max_num_atoms": max_num_atoms,
         "phelel_nosym": phelel_nosym,
+        "phonopy_max_num_atoms": phonopy_max_num_atoms,
         "phono3py_max_num_atoms": phono3py_max_num_atoms,
         "primitive_cell_choice": primitive_cell_choice,
         "symmetrize_cell": symmetrize_cell,
diff --git a/src/phelel/velph/cli/init/init.py b/src/phelel/velph/cli/init/init.py
index 51bad8c..85c24f5 100644
--- a/src/phelel/velph/cli/init/init.py
+++ b/src/phelel/velph/cli/init/init.py
@@ -614,35 +614,52 @@ def _get_toml_lines(
     phelel_dir_name: str = "phelel",
 ) -> list[str]:
     """Return velph-toml lines."""
+    # Parse [phelel] section for supercell dimension
     supercell_dimension = _get_supercell_dimension(
         velph_dict.get("phelel", {}),
         vip.max_num_atoms,
         sym_dataset,
         vip.find_primitive,
     )
-    if supercell_dimension is None:
-        click.echo("", err=True)
-        click.echo("Error | Supercell size could not be determined.", err=True)
-        click.echo(
-            "      | Specify max_num_atoms or [phelel.supercell_dimension].",
-            err=True,
-        )
-        return None
+    if supercell_dimension is not None:
+        click.echo("[phelel]")
+        _show_supercell_dimension(supercell_dimension)
+
+    # Parse [phonopy] section for supercell dimension
+    phonopy_supercell_dimension = _get_supercell_dimension(
+        velph_dict.get("phonopy", {}),
+        vip.phonopy_max_num_atoms,
+        sym_dataset,
+        vip.find_primitive,
+    )
+    if phonopy_supercell_dimension is not None:
+        click.echo("[phonopy]")
+        _show_supercell_dimension(phonopy_supercell_dimension)
 
+    # Parse [phono3py] section for supercell dimension
     phono3py_supercell_dimension = _get_supercell_dimension(
         velph_dict.get("phono3py", {}),
         vip.phono3py_max_num_atoms,
         sym_dataset,
         vip.find_primitive,
     )
-
-    click.echo("[phelel]")
-    _show_supercell_dimension(supercell_dimension)
-
     if phono3py_supercell_dimension is not None:
         click.echo("[phono3py]")
         _show_supercell_dimension(phono3py_supercell_dimension)
 
+    if (
+        supercell_dimension is None
+        and phonopy_supercell_dimension is None
+        and phono3py_supercell_dimension is None
+    ):
+        click.echo("", err=True)
+        click.echo("Error | Supercell size could not be determined.", err=True)
+        click.echo(
+            "      | Specify --max_num_atoms or [phelel.supercell_dimension].",
+            err=True,
+        )
+        return None
+
     (
         kpoints_dict,
         kpoints_dense_dict,
@@ -657,6 +674,7 @@ def _get_toml_lines(
         primitive,
         sym_dataset,
         supercell_dimension,
+        phonopy_supercell_dimension,
         phono3py_supercell_dimension,
         cell_choices["nac"],
         cell_choices["relax"],
@@ -696,6 +714,15 @@ def _get_toml_lines(
             vip.phelel_nosym,
         )
 
+    # [phonopy]
+    if phonopy_supercell_dimension is not None:
+        lines += ["[phonopy]"]
+        lines += _get_supercell_dimension_lines(phonopy_supercell_dimension)
+        lines += _get_displacement_settings_lines(
+            velph_dict, "phonopy", vip.amplitude, vip.diagonal, vip.plusminus
+        )
+        lines.append("")
+
     # [phono3py]
     if phono3py_supercell_dimension is not None:
         lines += ["[phono3py]"]
@@ -798,7 +825,8 @@ def _get_kpoints_dict(
     primitive: PhonopyAtoms,
     sym_dataset: SpglibDataset,
     supercell_dimension: np.ndarray,
-    phonon_supercell_dimension: Optional[np.ndarray],
+    phonopy_supercell_dimension: Optional[np.ndarray],
+    phono3py_supercell_dimension: Optional[np.ndarray],
     cell_for_nac: CellChoice,
     cell_for_relax: CellChoice,
     phelel_dir_name: str = "phelel",
@@ -824,29 +852,47 @@ def _get_kpoints_dict(
     )
 
     # Grid matrix for supercell
-    supercell = get_supercell(unitcell, supercell_dimension)
-    sym_dataset_super = get_symmetry_dataset(supercell, tolerance=vip_tolerance)
-    gm_super = GridMatrix(
-        2 * np.pi / vip_kspacing,
-        lattice=supercell.cell,
-        symmetry_dataset=sym_dataset_super,
-        use_grg=False,
-    )
+    if supercell_dimension is not None:
+        supercell = get_supercell(unitcell, supercell_dimension)
+        sym_dataset_super = get_symmetry_dataset(supercell, tolerance=vip_tolerance)
+        gm_super = GridMatrix(
+            2 * np.pi / vip_kspacing,
+            lattice=supercell.cell,
+            symmetry_dataset=sym_dataset_super,
+            use_grg=False,
+        )
+    else:
+        gm_super = None
+
+    # Grid matrix for phonopy supercell
+    if phonopy_supercell_dimension is not None:
+        phonopy_supercell = get_supercell(unitcell, phonopy_supercell_dimension)
+        sym_dataset_super = get_symmetry_dataset(
+            phonopy_supercell, tolerance=vip_tolerance
+        )
+        gm_phonopy_super = GridMatrix(
+            2 * np.pi / vip_kspacing,
+            lattice=phonopy_supercell.cell,
+            symmetry_dataset=sym_dataset_super,
+            use_grg=False,
+        )
+    else:
+        gm_phonopy_super = gm_super
 
-    # Grid matrix for phonon supercell
-    if phonon_supercell_dimension is not None:
-        phonon_supercell = get_supercell(unitcell, phonon_supercell_dimension)
+    # Grid matrix for phono3py supercell
+    if phono3py_supercell_dimension is not None:
+        phono3py_supercell = get_supercell(unitcell, phono3py_supercell_dimension)
         sym_dataset_super = get_symmetry_dataset(
-            phonon_supercell, tolerance=vip_tolerance
+            phono3py_supercell, tolerance=vip_tolerance
         )
-        gm_phonon_super = GridMatrix(
+        gm_phono3py_super = GridMatrix(
             2 * np.pi / vip_kspacing,
-            lattice=phonon_supercell.cell,
+            lattice=phono3py_supercell.cell,
             symmetry_dataset=sym_dataset_super,
             use_grg=False,
         )
     else:
-        gm_phonon_super = gm_super
+        gm_phono3py_super = gm_super
 
     # Dense grid matrix for primitive cell
     gm_dense_prim = GridMatrix(
@@ -861,12 +907,15 @@ def _get_kpoints_dict(
         gm,
         gm_prim,
         gm_super,
-        gm_phonon_super,
+        gm_phonopy_super,
+        gm_phono3py_super,
         cell_for_nac,
         cell_for_relax,
         phelel_dir_name=phelel_dir_name,
     )
-    kpoints_dense_dict = _get_kpoints_by_kspacing_dense(gm_dense_prim)
+    kpoints_dense_dict = _get_kpoints_by_kspacing_dense(
+        gm_dense_prim, with_phelel=phelel_dir_name in kpoints_dict
+    )
     qpoints_dict: dict = {}
     kpoints_opt_dict: dict = {}
 
@@ -952,7 +1001,8 @@ def _get_kpoints_by_kspacing(
     gm: GridMatrix,
     gm_prim: GridMatrix,
     gm_super: GridMatrix,
-    gm_phonon_super: Optional[GridMatrix],
+    gm_phonopy_super: Optional[GridMatrix],
+    gm_phono3py_super: Optional[GridMatrix],
     cell_for_nac: CellChoice,
     cell_for_relax: CellChoice,
     phelel_dir_name: str = "phelel",
@@ -966,7 +1016,11 @@ def _get_kpoints_by_kspacing(
     gm_prim : GridMatrix
         Grid matrix of primitive cell. The primitive cell can be a reduced cell.
     gm_super : GridMatrix
-        Grid matrix of supercell.
+        Grid matrix of phelel supercell.
+    gm_phonopy_super : GridMatrix
+        Grid matrix of phonopy supercell.
+    gm_phono3py_super : GridMatrix
+        Grid matrix of phono3py supercell.
     cell_for_nac : CellChoice
         Cell choice for NAC calculation among unit cell or primitive cell. The
         primitive cell can be a reduced cell.
@@ -981,17 +1035,27 @@ def _get_kpoints_by_kspacing(
         be calc_type names such as "phelel", "phelel.phonon", "selfenergy", etc.
 
     """
-    if gm_super.grid_matrix is None:
-        supercell_kpoints = {"mesh": gm_super.D_diag}
-    else:
-        supercell_kpoints = {"mesh": gm_super.grid_matrix}
-    if gm_phonon_super is not None:
+    if gm_super is not None:
         if gm_super.grid_matrix is None:
-            phonon_supercell_kpoints = {"mesh": gm_phonon_super.D_diag}
+            supercell_kpoints = {"mesh": gm_super.D_diag}
+        else:
+            supercell_kpoints = {"mesh": gm_super.grid_matrix}
+    else:
+        supercell_kpoints = None
+    if gm_phonopy_super is not None:
+        if gm_phonopy_super.grid_matrix is None:
+            phonopy_supercell_kpoints = {"mesh": gm_phonopy_super.D_diag}
+        else:
+            phonopy_supercell_kpoints = {"mesh": gm_phonopy_super.grid_matrix}
+    else:
+        phonopy_supercell_kpoints = None
+    if gm_phono3py_super is not None:
+        if gm_phono3py_super.grid_matrix is None:
+            phono3py_supercell_kpoints = {"mesh": gm_phono3py_super.D_diag}
         else:
-            phonon_supercell_kpoints = {"mesh": gm_phonon_super.grid_matrix}
+            phono3py_supercell_kpoints = {"mesh": gm_phono3py_super.grid_matrix}
     else:
-        phonon_supercell_kpoints = None
+        phono3py_supercell_kpoints = None
     if gm_prim.grid_matrix is None:
         selfenergy_kpoints = {"mesh": gm_prim.D_diag}
         el_bands_kpoints = {"mesh": gm_prim.D_diag}
@@ -1025,21 +1089,34 @@ def _get_kpoints_by_kspacing(
         relax_kpoints = {"mesh": gm_relax.grid_matrix}
 
     # keys are calc_types.
-    return {
+    kpoints_dict = {
         phelel_dir_name: supercell_kpoints,
-        "phono3py": supercell_kpoints,
-        "phelel.phonon": phonon_supercell_kpoints,
-        "phono3py.phonon": phonon_supercell_kpoints,
-        "selfenergy": selfenergy_kpoints,
-        "transport": selfenergy_kpoints,
+        "phonopy": phonopy_supercell_kpoints,
+        "phono3py": phono3py_supercell_kpoints,
         "relax": relax_kpoints,
         "nac": nac_kpoints,
         "el_bands": el_bands_kpoints,
         "ph_bands": ph_bands_kpoints,
     }
+    if gm_super:
+        kpoints_dict.update(
+            {
+                "selfenergy": selfenergy_kpoints,
+                "transport": selfenergy_kpoints,
+            }
+        )
+
+    return_kpoints_dict = {}
+    for key, val in kpoints_dict.items():
+        if val is not None:
+            return_kpoints_dict[key] = val
 
+    return return_kpoints_dict
 
-def _get_kpoints_by_kspacing_dense(gm_dense_prim: GridMatrix) -> dict:
+
+def _get_kpoints_by_kspacing_dense(
+    gm_dense_prim: GridMatrix, with_phelel: bool = True
+) -> dict:
     """Return kpoints dict of dense grid.
 
     Returns
@@ -1055,12 +1132,17 @@ def _get_kpoints_by_kspacing_dense(gm_dense_prim: GridMatrix) -> dict:
     else:
         selfenergy_kpoints_dense = {"mesh": gm_dense_prim.grid_matrix}
         el_bands_kpoints_dense = {"mesh": gm_dense_prim.grid_matrix}
-    # keys are calc_types.
-    return {
-        "selfenergy": selfenergy_kpoints_dense,
-        "transport": selfenergy_kpoints_dense,
-        "el_bands": el_bands_kpoints_dense,
-    }
+
+    if with_phelel:
+        return {
+            "selfenergy": selfenergy_kpoints_dense,
+            "transport": selfenergy_kpoints_dense,
+            "el_bands": el_bands_kpoints_dense,
+        }
+    else:
+        return {
+            "el_bands": el_bands_kpoints_dense,
+        }
 
 
 def _get_vasp_lines(
@@ -1077,8 +1159,8 @@ def _get_vasp_lines(
 
     lines = []
 
-    for calc_type in (phelel_dir_name, "phono3py"):
-        if calc_type in vasp_dict:
+    for calc_type in (phelel_dir_name, "phonopy", "phono3py"):
+        if calc_type in vasp_dict and calc_type in kpoints_dict:
             _vasp_dict = vasp_dict[calc_type]
             _add_incar_lines(lines, _vasp_dict, incar_commons, calc_type)
             lines.append(f"[vasp.{calc_type}.kpoints]")
@@ -1086,20 +1168,12 @@ def _get_vasp_lines(
             _add_calc_type_scheduler_lines(lines, _vasp_dict, calc_type)
             lines.append("")
 
-    for calc_type in (phelel_dir_name, "phono3py"):
-        if calc_type in vasp_dict:
-            if "phonon" in vasp_dict[calc_type]:
-                _vasp_dict = vasp_dict[calc_type]["phonon"]
-                _add_incar_lines(
-                    lines, _vasp_dict, incar_commons, f"{calc_type}.phonon"
-                )
-                lines.append(f"[vasp.{calc_type}.phonon.kpoints]")
-                _add_kpoints_lines(lines, kpoints_dict[f"{calc_type}.phonon"])
-                _add_calc_type_scheduler_lines(lines, _vasp_dict, f"{calc_type}.phonon")
-                lines.append("")
-
     for calc_type in ("selfenergy", "transport"):
-        if calc_type in vasp_dict:
+        if (
+            calc_type in vasp_dict
+            and calc_type in kpoints_dict
+            and calc_type in kpoints_dense_dict
+        ):
             # primitive cell
             _vasp_dict = vasp_dict[calc_type]
             _add_incar_lines(lines, _vasp_dict, incar_commons, calc_type)
@@ -1202,22 +1276,24 @@ def _get_phelel_lines(
     lines.append("[phelel]")
     lines.append(f'version = "{__version__}"')
 
-    lines += _get_supercell_dimension_lines(supercell_dimension)
-    lines += _get_displacement_settings_lines(
-        velph_dict, "phelel", amplitude, diagonal, plusminus
-    )
+    if supercell_dimension is not None:
+        lines += _get_supercell_dimension_lines(supercell_dimension)
+        lines += _get_displacement_settings_lines(
+            velph_dict, "phelel", amplitude, diagonal, plusminus
+        )
 
-    if phelel_nosym:
-        lines.append("nosym = true")
+        if phelel_nosym:
+            lines.append("nosym = true")
+
+        fft_mesh = _get_fft_mesh(velph_dict, primitive)
+        try:
+            if fft_mesh is None and "fft_mesh" in velph_dict["phelel"]:
+                fft_mesh = velph_dict["phelel"]["fft_mesh"]
+            if fft_mesh is not None:
+                lines.append("fft_mesh = [{:d}, {:d}, {:d}]".format(*fft_mesh))
+        except KeyError:
+            pass
 
-    fft_mesh = _get_fft_mesh(velph_dict, primitive)
-    try:
-        if fft_mesh is None and "fft_mesh" in velph_dict["phelel"]:
-            fft_mesh = velph_dict["phelel"]["fft_mesh"]
-        if fft_mesh is not None:
-            lines.append("fft_mesh = [{:d}, {:d}, {:d}]".format(*fft_mesh))
-    except KeyError:
-        pass
     lines.append("")
     return lines
 
@@ -1333,7 +1409,7 @@ def _get_supercell_dimension_lines(
 
 def _get_displacement_settings_lines(
     velph_dict: dict,
-    calc_type: Literal["phelel", "phono3py"],
+    calc_type: Literal["phelel", "phonopy", "phono3py"],
     amplitude: Optional[float],
     diagonal: Optional[bool],
     plusminus: Optional[bool],
diff --git a/src/phelel/velph/cli/phono3py/__init__.py b/src/phelel/velph/cli/phono3py/__init__.py
index 446bfeb..c9006dc 100644
--- a/src/phelel/velph/cli/phono3py/__init__.py
+++ b/src/phelel/velph/cli/phono3py/__init__.py
@@ -1,4 +1,4 @@
-"""velph command line tool / velph-supercell."""
+"""velph command line tool / velph-phono3py."""
 
 from __future__ import annotations
 
@@ -38,36 +38,27 @@ def cmd_phono3py():
     type=int,
     help="Number of snapshots of supercells with random directional displacement.",
 )
-@click.option(
-    "--rd-fc2",
-    "random_displacements_fc2",
-    nargs=1,
-    default=None,
-    type=int,
-    help=(
-        "Number of snapshots of phonon supercells "
-        "with random directional displacement."
-    ),
-)
 @click.help_option("-h", "--help")
 def cmd_init(
     toml_filename: str,
     random_displacements: Optional[int],
-    random_displacements_fc2: Optional[int],
 ):
-    """Generate displacements and write phelel_disp.yaml."""
+    """Generate displacements and write phono3py_disp.yaml."""
     with open(toml_filename, "rb") as f:
         toml_dict = tomli.load(f)
 
-    ph3py = run_init(
+    if "phono3py" not in toml_dict:
+        click.echo(f'[phono3py] section not found in "{toml_filename}" file.', err=True)
+        return
+
+    ph3 = run_init(
         toml_dict,
         number_of_snapshots=random_displacements,
-        number_of_snapshots_fc2=random_displacements_fc2,
     )
 
     phono3py_yaml_filename = pathlib.Path("phono3py/phono3py_disp.yaml")
     phono3py_yaml_filename.parent.mkdir(parents=True, exist_ok=True)
-    ph3py.save(phono3py_yaml_filename)
+    ph3.save(phono3py_yaml_filename)
 
     click.echo(f'"{phono3py_yaml_filename}" was generated. ')
     click.echo('VASP input files will be generated by "velph phono3py generate".')
diff --git a/src/phelel/velph/cli/phono3py/init.py b/src/phelel/velph/cli/phono3py/init.py
index 41be45c..455e25b 100644
--- a/src/phelel/velph/cli/phono3py/init.py
+++ b/src/phelel/velph/cli/phono3py/init.py
@@ -18,7 +18,6 @@ def run_init(
     toml_dict: dict,
     current_directory: pathlib.Path = pathlib.Path(""),
     number_of_snapshots: Optional[int] = None,
-    number_of_snapshots_fc2: Optional[int] = None,
 ) -> Optional[Phono3py]:
     """Generate displacements and write phono3py_disp.yaml.
 
@@ -26,13 +25,11 @@ def run_init(
         Used for test.
 
     """
+    if "phono3py" not in toml_dict:
+        raise RuntimeError("[phono3py] section not found in toml file.")
+
     convcell = parse_cell_dict(toml_dict["unitcell"])
-    supercell_matrix = toml_dict["phelel"].get("supercell_dimension", None)
-    if "phono3py" in toml_dict:
-        phonon_supercell_matrix = supercell_matrix
-        supercell_matrix = toml_dict["phono3py"].get("supercell_dimension", None)
-    else:
-        phonon_supercell_matrix = None
+    supercell_matrix = toml_dict["phono3py"].get("supercell_dimension", None)
     if "primitive_cell" in toml_dict:
         primitive = parse_cell_dict(toml_dict["primitive_cell"])
         primitive_matrix = np.dot(np.linalg.inv(convcell.cell.T), primitive.cell.T)
@@ -42,7 +39,7 @@ def run_init(
 
     is_symmetry = True
     try:
-        if toml_dict["phelel"]["nosym"] is True:
+        if toml_dict["phono3py"]["nosym"] is True:
             is_symmetry = False
     except KeyError:
         pass
@@ -50,17 +47,15 @@ def run_init(
     ph3py = Phono3py(
         convcell,
         supercell_matrix=supercell_matrix,
-        phonon_supercell_matrix=phonon_supercell_matrix,
         primitive_matrix=primitive_matrix,
         is_symmetry=is_symmetry,
         calculator="vasp",
     )
 
-    amplitude = toml_dict["phelel"].get("amplitude", None)
-
+    amplitude = toml_dict["phono3py"].get("amplitude", None)
     if number_of_snapshots is None:
-        is_diagonal = toml_dict["phelel"].get("diagonal", True)
-        is_plusminus = toml_dict["phelel"].get("plusminus", "auto")
+        is_diagonal = toml_dict["phono3py"].get("diagonal", True)
+        is_plusminus = toml_dict["phono3py"].get("plusminus", "auto")
     else:
         is_diagonal = False
         is_plusminus = False
@@ -72,7 +67,6 @@ def run_init(
         is_plusminus=is_plusminus,
         is_diagonal=is_diagonal,
         number_of_snapshots=number_of_snapshots,
-        number_of_snapshots_fc2=number_of_snapshots_fc2,
     )
 
     nac_directory = current_directory / "nac"
diff --git a/src/phelel/velph/cli/utils.py b/src/phelel/velph/cli/utils.py
index 3f83c0a..837e711 100644
--- a/src/phelel/velph/cli/utils.py
+++ b/src/phelel/velph/cli/utils.py
@@ -84,6 +84,7 @@ class VelphInitParams:
     kspacing_dense: Optional[float] = 0.05
     magmom: Optional[str] = None
     max_num_atoms: Optional[int] = None
+    phonopy_max_num_atoms: Optional[int] = None
     phono3py_max_num_atoms: Optional[int] = None
     phelel_nosym: Optional[bool] = False
     primitive_cell_choice: Optional[PrimitiveCellChoice] = (
diff --git a/test/velph/cli/phono3py/init/test_phono3py_init.py b/test/velph/cli/phono3py/init/test_phono3py_init.py
index b8c3b8f..ef9866e 100644
--- a/test/velph/cli/phono3py/init/test_phono3py_init.py
+++ b/test/velph/cli/phono3py/init/test_phono3py_init.py
@@ -39,23 +39,13 @@ def test_phono3py_init_random_displacements(distance: float):
 coordinates = [  0.666666666666664,  0.333333333333336,  0.750000000000000 ]
 magnetic_moment = [ 0.00000000, 0.00000000, 0.00000000 ]"""
 
-    toml_dict = tomli.loads(phelel_str + unitcell_str)
-    ph3 = run_init(toml_dict, number_of_snapshots=10)
-    np.testing.assert_array_equal(ph3.supercell_matrix, np.diag([4, 4, 2]))
-    assert len(ph3.supercell) == 64
-    assert ph3.displacements.shape == (10, 64, 3)
-    np.testing.assert_allclose(np.linalg.norm(ph3.displacements, axis=2), distance)
-
-    toml_dict = tomli.loads(phelel_str + phono3py_str + unitcell_str)
-    ph3 = run_init(toml_dict, number_of_snapshots=10, number_of_snapshots_fc2=4)
-    np.testing.assert_array_equal(ph3.supercell_matrix, np.diag([2, 2, 1]))
-    np.testing.assert_array_equal(ph3.phonon_supercell_matrix, np.diag([4, 4, 2]))
-
-    assert len(ph3.supercell) == 8
-    assert ph3.displacements.shape == (10, 8, 3)
-    assert len(ph3.phonon_supercell) == 64
-    assert ph3.phonon_displacements.shape == (4, 64, 3)
-    np.testing.assert_allclose(np.linalg.norm(ph3.displacements, axis=2), distance)
-    np.testing.assert_allclose(
-        np.linalg.norm(ph3.phonon_displacements, axis=2), distance
-    )
+    for toml_dict in (
+        tomli.loads(phono3py_str + unitcell_str),
+        tomli.loads(phelel_str + phono3py_str + unitcell_str),
+    ):
+        ph3 = run_init(toml_dict, number_of_snapshots=10)
+        np.testing.assert_array_equal(ph3.supercell_matrix, np.diag([2, 2, 1]))
+
+        assert len(ph3.supercell) == 8
+        assert ph3.displacements.shape == (10, 8, 3)
+        np.testing.assert_allclose(np.linalg.norm(ph3.displacements, axis=2), distance)

From 7e96636a511685a7735b86aa4e7788228350d618 Mon Sep 17 00:00:00 2001
From: Atsushi Togo <atz.togo@gmail.com>
Date: Mon, 9 Dec 2024 15:31:44 +0900
Subject: [PATCH 2/4] Add velph-phonopy files

---
 src/phelel/velph/cli/phonopy/__init__.py |  86 +++++++++++++++++
 src/phelel/velph/cli/phonopy/generate.py |  36 ++++++++
 src/phelel/velph/cli/phonopy/init.py     | 113 +++++++++++++++++++++++
 3 files changed, 235 insertions(+)
 create mode 100644 src/phelel/velph/cli/phonopy/__init__.py
 create mode 100644 src/phelel/velph/cli/phonopy/generate.py
 create mode 100644 src/phelel/velph/cli/phonopy/init.py

diff --git a/src/phelel/velph/cli/phonopy/__init__.py b/src/phelel/velph/cli/phonopy/__init__.py
new file mode 100644
index 0000000..8deab3e
--- /dev/null
+++ b/src/phelel/velph/cli/phonopy/__init__.py
@@ -0,0 +1,86 @@
+"""velph command line tool / velph-phonopy."""
+
+from __future__ import annotations
+
+import pathlib
+from typing import Optional
+
+import click
+import tomli
+
+from phelel.velph.cli import cmd_root
+from phelel.velph.cli.phonopy.generate import write_supercell_input_files
+from phelel.velph.cli.phonopy.init import run_init
+
+
+@cmd_root.group("phonopy")
+@click.help_option("-h", "--help")
+def cmd_phonopy():
+    """Choose phonopy options."""
+    pass
+
+
+#
+# velph phonopy init
+#
+@cmd_phonopy.command("init")
+@click.argument(
+    "toml_filename",
+    nargs=1,
+    type=click.Path(),
+    default="velph.toml",
+)
+@click.option(
+    "--rd",
+    "random_displacements",
+    nargs=1,
+    default=None,
+    type=int,
+    help="Number of snapshots of supercells with random directional displacement.",
+)
+@click.help_option("-h", "--help")
+def cmd_init(
+    toml_filename: str,
+    random_displacements: Optional[int],
+):
+    """Generate displacements and write phonopy_disp.yaml."""
+    with open(toml_filename, "rb") as f:
+        toml_dict = tomli.load(f)
+
+    if "phonopy" not in toml_dict:
+        click.echo(f'[phonopy] section not found in "{toml_filename}" file.', err=True)
+        return
+
+    ph = run_init(
+        toml_dict,
+        number_of_snapshots=random_displacements,
+    )
+
+    phonopy_yaml_filename = pathlib.Path("phonopy/phonopy_disp.yaml")
+    phonopy_yaml_filename.parent.mkdir(parents=True, exist_ok=True)
+    ph.save(phonopy_yaml_filename)
+
+    click.echo(f'"{phonopy_yaml_filename}" was generated. ')
+    click.echo('VASP input files will be generated by "velph phonopy generate".')
+
+
+#
+# velph phonopy generate
+#
+@cmd_phonopy.command("generate")
+@click.argument(
+    "toml_filename",
+    nargs=1,
+    type=click.Path(),
+    default="velph.toml",
+)
+@click.help_option("-h", "--help")
+def cmd_generate(toml_filename: str):
+    """Generate phonopy input files."""
+    yaml_filename = "phonopy/phonopy_disp.yaml"
+    if not pathlib.Path("POTCAR").exists():
+        click.echo('"POTCAR" not found in current directory.')
+
+    write_supercell_input_files(
+        pathlib.Path(toml_filename), pathlib.Path(yaml_filename)
+    )
diff --git a/src/phelel/velph/cli/phonopy/generate.py b/src/phelel/velph/cli/phonopy/generate.py
new file mode 100644
index 0000000..c3e9dea
--- /dev/null
+++ b/src/phelel/velph/cli/phonopy/generate.py
@@ -0,0 +1,36 @@
+"""Implementation of velph-phonopy-generate."""
+
+from __future__ import annotations
+
+import pathlib
+
+import click
+import phono3py
+import tomli
+
+from phelel.velph.cli.phelel.generate import (
+    write_phonon_supercells,
+    write_supercells,
+)
+
+
+def write_supercell_input_files(
+    toml_filename: pathlib.Path,
+    phono3py_yaml_filename: pathlib.Path,
+) -> None:
+    """Generate supercells."""
+    if not phono3py_yaml_filename.exists():
+        click.echo(f'File "{phono3py_yaml_filename}" not found.', err=True)
+        click.echo('Run "velph phono3py init" if necessary.', err=True)
+        return None
+
+    ph3py = phono3py.load(phono3py_yaml_filename)
+    with open(toml_filename, "rb") as f:
+        toml_dict = tomli.load(f)
+
+    write_supercells(ph3py, toml_dict, dir_name="phono3py")
+    if ph3py.phonon_supercell_matrix is not None:
+        if "phonon" in toml_dict["vasp"]["phono3py"]:
+            write_phonon_supercells(ph3py, toml_dict, dir_name="phono3py")
+        else:
+            print(f'[vasp.phono3py.phonon.*] not found in "{toml_filename}"')
diff --git a/src/phelel/velph/cli/phonopy/init.py b/src/phelel/velph/cli/phonopy/init.py
new file mode 100644
index 0000000..62fe76f
--- /dev/null
+++ b/src/phelel/velph/cli/phonopy/init.py
@@ -0,0 +1,113 @@
+"""Implementation of velph-phonopy-init."""
+
+from __future__ import annotations
+
+import pathlib
+from typing import Optional, Union
+
+import click
+import numpy as np
+from phonopy import Phonopy
+from phonopy.interface.calculator import get_default_displacement_distance
+from phonopy.structure.atoms import parse_cell_dict
+
+from phelel.velph.cli.utils import get_nac_params
+
+
+def run_init(
+    toml_dict: dict,
+    current_directory: pathlib.Path = pathlib.Path(""),
+    number_of_snapshots: Optional[int] = None,
+) -> Optional[Phonopy]:
+    """Generate displacements and write phonopy_disp.yaml.
+
+    current_directory : Path
+        Used for test.
+
+    """
+    if "phonopy" not in toml_dict:
+        raise RuntimeError("[phonopy] section not found in toml file.")
+
+    convcell = parse_cell_dict(toml_dict["unitcell"])
+    supercell_matrix = toml_dict["phonopy"].get("supercell_dimension", None)
+    if "primitive_cell" in toml_dict:
+        primitive = parse_cell_dict(toml_dict["primitive_cell"])
+        primitive_matrix = np.dot(np.linalg.inv(convcell.cell.T), primitive.cell.T)
+    else:
+        primitive = convcell
+        primitive_matrix = None
+
+    is_symmetry = True
+    try:
+        if toml_dict["phonopy"]["nosym"] is True:
+            is_symmetry = False
+    except KeyError:
+        pass
+
+    ph = Phonopy(
+        convcell,
+        supercell_matrix=supercell_matrix,
+        primitive_matrix=primitive_matrix,
+        is_symmetry=is_symmetry,
+        calculator="vasp",
+    )
+
+    amplitude = toml_dict["phonopy"].get("amplitude", None)
+    if number_of_snapshots is None:
+        is_diagonal = toml_dict["phonopy"].get("diagonal", True)
+        is_plusminus = toml_dict["phonopy"].get("plusminus", "auto")
+    else:
+        is_diagonal = False
+        is_plusminus = False
+
+    _generate_phonopy_supercells(
+        ph,
+        interface_mode="vasp",
+        distance=amplitude,
+        is_plusminus=is_plusminus,
+        is_diagonal=is_diagonal,
+        number_of_snapshots=number_of_snapshots,
+    )
+
+    nac_directory = current_directory / "nac"
+    if nac_directory.exists():
+        click.echo('Found "nac" directory. Read NAC params.')
+        vasprun_path = nac_directory / "vasprun.xml"
+        if vasprun_path.exists():
+            nac_params = get_nac_params(
+                toml_dict,
+                vasprun_path,
+                primitive,
+                convcell,
+                is_symmetry,
+            )
+            if nac_params is not None:
+                ph.nac_params = nac_params
+        else:
+            click.echo('Not found "nac/vasprun.xml". NAC params were not included.')
+
+    return ph
+
+
+def _generate_phonopy_supercells(
+    phonopy: Phonopy,
+    interface_mode: str = "vasp",
+    distance: Optional[float] = None,
+    is_plusminus: Union[str, bool] = "auto",
+    is_diagonal: bool = True,
+    number_of_snapshots: Optional[int] = None,
+):
+    """Generate phelel supercells."""
+    if distance is None:
+        _distance = get_default_displacement_distance(interface_mode)
+    else:
+        _distance = distance
+
+    phonopy.generate_displacements(
+        distance=_distance,
+        is_plusminus=is_plusminus,
+        is_diagonal=is_diagonal,
+        number_of_snapshots=number_of_snapshots,
+    )
+    click.echo(f"Displacement distance: {_distance}")
+    click.echo(f"Number of displacements: {len(phonopy.supercells_with_displacements)}")

From 28181ac75eb4379cd6417abc3fc1c32f856a93ab Mon Sep 17 00:00:00 2001
From: Atsushi Togo <atz.togo@gmail.com>
Date: Mon, 9 Dec 2024 15:37:07 +0900
Subject: [PATCH 3/4] Fix velph-phonopy-generate

---
 src/phelel/velph/cli/phelel/generate.py  |  3 ++-
 src/phelel/velph/cli/phonopy/generate.py | 20 +++++++-------------
 2 files changed, 9 insertions(+), 14 deletions(-)

diff --git a/src/phelel/velph/cli/phelel/generate.py b/src/phelel/velph/cli/phelel/generate.py
index d91dbad..945f0b5 100644
--- a/src/phelel/velph/cli/phelel/generate.py
+++ b/src/phelel/velph/cli/phelel/generate.py
@@ -9,6 +9,7 @@
 import click
 import tomli
 from phono3py import Phono3py
+from phonopy import Phonopy
 from phonopy.interface.calculator import write_crystal_structure
 
 import phelel
@@ -47,7 +48,7 @@ def write_supercell_input_files(
 
 
 def write_supercells(
-    phe: Union[Phelel, Phono3py], toml_dict: dict, dir_name: str = "phelel"
+    phe: Union[Phelel, Phonopy, Phono3py], toml_dict: dict, dir_name: str = "phelel"
 ):
     """Write VASP input for supercells.
 
diff --git a/src/phelel/velph/cli/phonopy/generate.py b/src/phelel/velph/cli/phonopy/generate.py
index c3e9dea..6fc89aa 100644
--- a/src/phelel/velph/cli/phonopy/generate.py
+++ b/src/phelel/velph/cli/phonopy/generate.py
@@ -5,32 +5,26 @@
 import pathlib
 
 import click
-import phono3py
+import phonopy
 import tomli
 
 from phelel.velph.cli.phelel.generate import (
-    write_phonon_supercells,
     write_supercells,
 )
 
 
 def write_supercell_input_files(
     toml_filename: pathlib.Path,
-    phono3py_yaml_filename: pathlib.Path,
+    phonopy_yaml_filename: pathlib.Path,
 ) -> None:
     """Generate supercells."""
-    if not phono3py_yaml_filename.exists():
-        click.echo(f'File "{phono3py_yaml_filename}" not found.', err=True)
-        click.echo('Run "velph phono3py init" if necessary.', err=True)
+    if not phonopy_yaml_filename.exists():
+        click.echo(f'File "{phonopy_yaml_filename}" not found.', err=True)
+        click.echo('Run "velph phonopy init" if necessary.', err=True)
         return None
 
-    ph3py = phono3py.load(phono3py_yaml_filename)
+    ph = phonopy.load(phonopy_yaml_filename)
     with open(toml_filename, "rb") as f:
         toml_dict = tomli.load(f)
 
-    write_supercells(ph3py, toml_dict, dir_name="phono3py")
-    if ph3py.phonon_supercell_matrix is not None:
-        if "phonon" in toml_dict["vasp"]["phono3py"]:
-            write_phonon_supercells(ph3py, toml_dict, dir_name="phono3py")
-        else:
-            print(f'[vasp.phono3py.phonon.*] not found in "{toml_filename}"')
+    write_supercells(ph, toml_dict, dir_name="phonopy")

From 007daed7226b33d9cf5f9d1fd87e4ffc36b70af0 Mon Sep 17 00:00:00 2001
From: Atsushi Togo <atz.togo@gmail.com>
Date: Mon, 9 Dec 2024 15:41:05 +0900
Subject: [PATCH 4/4] Add test for velph-phonopy-init

---
 test/velph/cli/__init__.py                    |  1 -
 .../cli/phonopy/init/test_phonopy_init.py     | 51 +++++++++++++++++++
 2 files changed, 51 insertions(+), 1 deletion(-)
 delete mode 100644 test/velph/cli/__init__.py
 create mode 100644 test/velph/cli/phonopy/init/test_phonopy_init.py

diff --git a/test/velph/cli/__init__.py b/test/velph/cli/__init__.py
deleted file mode 100644
index 31c1276..0000000
--- a/test/velph/cli/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-"""CLI implementation."""
diff --git a/test/velph/cli/phonopy/init/test_phonopy_init.py b/test/velph/cli/phonopy/init/test_phonopy_init.py
new file mode 100644
index 0000000..34bb628
--- /dev/null
+++ b/test/velph/cli/phonopy/init/test_phonopy_init.py
@@ -0,0 +1,51 @@
+"""Tests velph-phonopy-init."""
+
+import numpy as np
+import pytest
+import tomli
+
+from phelel.velph.cli.phonopy.init import run_init
+
+
+@pytest.mark.parametrize("distance", [0.03, 0.05])
+def test_phonopy_init_random_displacements(distance: float):
+    """Test of plusminus and diagonal with Ti."""
+    phelel_str = f"""title = "VASP el-ph settings"
+
+[phelel]
+supercell_dimension = [4, 4, 2]
+amplitude = {distance}
+fft_mesh = [18, 18, 28]
+"""
+
+    phonopy_str = f"""[phonopy]
+supercell_dimension = [2, 2, 1]
+amplitude = {distance}
+"""
+
+    unitcell_str = """
+[unitcell]
+lattice = [
+  [     2.930720886111760,     0.000000000000000,     0.000000000000000 ], # a
+  [    -1.465360443055880,     2.538078738774425,     0.000000000000000 ], # b
+  [     0.000000000000000,     0.000000000000000,     4.646120482318025 ], # c
+]
+[[unitcell.points]]  # 1
+symbol = "Ti"
+coordinates = [  0.333333333333336,  0.666666666666664,  0.250000000000000 ]
+magnetic_moment = [ 0.00000000, 0.00000000, 0.00000000 ]
+[[unitcell.points]]  # 2
+symbol = "Ti"
+coordinates = [  0.666666666666664,  0.333333333333336,  0.750000000000000 ]
+magnetic_moment = [ 0.00000000, 0.00000000, 0.00000000 ]"""
+
+    for toml_dict in (
+        tomli.loads(phonopy_str + unitcell_str),
+        tomli.loads(phelel_str + phonopy_str + unitcell_str),
+    ):
+        ph = run_init(toml_dict, number_of_snapshots=10)
+        np.testing.assert_array_equal(ph.supercell_matrix, np.diag([2, 2, 1]))
+
+        assert len(ph.supercell) == 8
+        assert ph.displacements.shape == (10, 8, 3)
+        np.testing.assert_allclose(np.linalg.norm(ph.displacements, axis=2), distance)