Include velph in phelel

.github/workflows/tests.yml

@@ -26,7 +26,7 @@ jobs:
       run: |
         conda activate test
         conda install --yes -c conda-forge python=${{ matrix.python-version }}
-        conda install --yes -c conda-forge phonopy phono3py finufft pytest codecov pytest-cov
+        conda install --yes -c conda-forge phonopy phono3py finufft tomli tomli-w seekpath pytest click codecov pytest-cov
     - name: Setup phelel
       run: |
         conda activate test

pyproject.toml

@@ -13,11 +13,16 @@ authors = [
 requires-python = ">=3.9"
 dependencies = [
     "phono3py",
+    "finufft",
+    "click",
+    "tomli",
+    "tomli-w"
 ]
 license = {file = "LICENSE"}
 
 [project.scripts]
 phelel = "phelel.scripts.phelel:run"
+velph = "velph.cli:cmd_root"
 
 [tool.setuptools.dynamic]
 version = { attr = "phelel.version.__version__" }
@@ -35,11 +40,6 @@ lint.extend-ignore = [
     "D417",
     "D100",
 ]
-exclude = [
-  "test/phonon/test_irreps.py",
-  "test/qha/test_electron.py",
-  "phonopy/interface/cp2k.py",
-]
 
 [tool.ruff.lint.pydocstyle]
 convention = "numpy"

src/phelel/velph/__init__.py

@@ -0,0 +1,5 @@
+"""velph module."""
+
+from phelel.velph import cli  # noqa F401
+from phelel.velph import templates  # noqa F401
+from phelel.velph import utils  # noqa F401

src/phelel/velph/cli/__init__.py

@@ -0,0 +1,24 @@
+"""velph command line tool module."""
+
+import click
+
+
+@click.group()
+@click.help_option("-h", "--help")
+def cmd_root():
+    """Command-line utility to help VASP el-ph calculation."""
+    pass
+
+
+from phelel.velph.cli.el_bands import cmd_el_bands  # noqa F401
+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.nac import cmd_nac  # noqa F401
+from phelel.velph.cli.ph_bands import cmd_ph_bands  # noqa F401
+from phelel.velph.cli.relax import cmd_relax  # noqa F401
+
+# The followings are written here to avoid ciclic import but needed.
+from phelel.velph.cli.selfenergy import cmd_selfenergy  # noqa F401
+from phelel.velph.cli.supercell import cmd_supercell  # noqa F401
+from phelel.velph.cli.transport import cmd_transport  # noqa F401

src/phelel/velph/cli/el_bands/__init__.py

@@ -0,0 +1,66 @@
+"""velph command line tool / velph-el_bands."""
+
+from __future__ import annotations
+
+import pathlib
+
+import click
+
+from phelel.velph.cli import cmd_root
+from phelel.velph.cli.el_bands.generate import write_input_files
+from phelel.velph.cli.el_bands.plot import plot_el_bandstructures
+
+
+@cmd_root.group("el_bands")
+@click.help_option("-h", "--help")
+def cmd_el_bands():
+    """Choose electronic band structure options."""
+    pass
+
+
+#
+# velph el_bands generate
+#
+@cmd_el_bands.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 input files to plot electronic band structure."""
+    if not pathlib.Path("POTCAR").exists():
+        click.echo('"POTCAR" not found in current directory.')
+
+    write_input_files(pathlib.Path(toml_filename))
+
+
+#
+# velph el_bands plot
+#
+@cmd_el_bands.command("plot")
+@click.option(
+    "--window",
+    required=True,
+    type=(float, float),
+    help="Energy window, emin and emax with respect to Fermi level.",
+)
+@click.help_option("-h", "--help")
+def cmd_plot(window: tuple[float, float]):
+    """Plot electronic band structure."""
+    vaspout_filename_bands = pathlib.Path("el_bands/bands/vaspout.h5")
+    if vaspout_filename_bands.exists():
+        vaspout_filename_dos = pathlib.Path("el_bands/dos/vaspout.h5")
+        if vaspout_filename_bands.exists():
+            click.echo(f'Found "{vaspout_filename_dos}". DOS will be plotted.')
+            plot_el_bandstructures(
+                window,
+                vaspout_filename_bands,
+                vaspout_filename_dos=vaspout_filename_dos,
+            )
+        else:
+            plot_el_bandstructures(window, vaspout_filename_bands)
+    else:
+        click.echo(f'"{vaspout_filename_bands}" not found.')

src/phelel/velph/cli/el_bands/generate.py

@@ -0,0 +1,101 @@
+"""Implementation of velph-relax-generate."""
+
+import copy
+import pathlib
+import shutil
+
+import click
+import tomli
+from phonopy.interface.calculator import write_crystal_structure
+from phonopy.structure.atoms import parse_cell_dict
+
+from phelel.velph.cli.utils import (
+    assert_kpoints_mesh_symmetry,
+    get_scheduler_dict,
+    write_incar,
+    write_kpoints_line_mode,
+    write_kpoints_mesh_mode,
+    write_launch_script,
+)
+
+
+def write_input_files(toml_filename: pathlib.Path) -> None:
+    """Generate VASP inputs to generate electronc band structure."""
+    with open(toml_filename, "rb") as f:
+        toml_dict = tomli.load(f)
+
+    main_directory_name = "el_bands"
+    main_directory = pathlib.Path(main_directory_name)
+    main_directory.mkdir(parents=True, exist_ok=True)
+
+    for calc_type in ("bands", "dos"):
+        directory_name = f"{main_directory_name}/{calc_type}"
+        directory = pathlib.Path(directory_name)
+        directory.mkdir(parents=True, exist_ok=True)
+
+        # POSCAR
+        primitive = parse_cell_dict(toml_dict["primitive_cell"])
+        write_crystal_structure(directory / "POSCAR", primitive)
+
+        # INCAR
+        if calc_type == "bands":
+            incar_dict = copy.deepcopy(toml_dict["vasp"]["el_bands"]["incar"])
+            for key in toml_dict["vasp"]["el_bands"]["incar"]:
+                if key.lower().strip() in ("lorbit", "nedos"):
+                    del incar_dict[key]
+            write_incar(incar_dict, directory, cell=primitive)
+        elif calc_type == "dos":
+            write_incar(
+                toml_dict["vasp"]["el_bands"]["incar"], directory, cell=primitive
+            )
+
+        # KPOINTS
+        kpoints_dict = toml_dict["vasp"]["el_bands"]["kpoints"]
+        assert_kpoints_mesh_symmetry(toml_dict, kpoints_dict, primitive)
+        write_kpoints_mesh_mode(
+            toml_dict["vasp"]["el_bands"]["incar"],
+            directory,
+            "vasp.el_bands.kpoints",
+            toml_dict["vasp"]["el_bands"]["kpoints"],
+        )
+
+        # KPOINTS_OPT
+        if calc_type == "bands":
+            if "path" in toml_dict["vasp"]["el_bands"]["kpoints_opt"]:
+                click.echo(
+                    "Seek-path (https://github.com/giovannipizzi/seekpath) is used."
+                )
+
+            write_kpoints_line_mode(
+                primitive,
+                directory,
+                "vasp.el_bands.kpoints_opt",
+                toml_dict["vasp"]["el_bands"]["kpoints_opt"],
+                kpoints_filename="KPOINTS_OPT",
+            )
+        elif calc_type == "dos":
+            if "kpoints_dense" not in toml_dict["vasp"]["el_bands"]:
+                raise RuntimeError(
+                    "[vasp.el_bands.kpoints_dense] section is necessary "
+                    "for electronic DOS calculation."
+                )
+            kpoints_dense_dict = toml_dict["vasp"]["el_bands"]["kpoints_dense"]
+            assert_kpoints_mesh_symmetry(toml_dict, kpoints_dense_dict, primitive)
+            write_kpoints_mesh_mode(
+                toml_dict["vasp"]["el_bands"]["incar"],
+                directory,
+                "vasp.el_bands.kpoints",
+                kpoints_dense_dict,
+                kpoints_filename="KPOINTS_OPT",
+            )
+
+        # POTCAR
+        if pathlib.Path("POTCAR").exists():
+            shutil.copy2(pathlib.Path("POTCAR"), directory / "POTCAR")
+
+        # Scheduler launch script
+        if "scheduler" in toml_dict:
+            scheduler_dict = get_scheduler_dict(toml_dict, "el_bands")
+            write_launch_script(scheduler_dict, directory)
+
+        click.echo(f'VASP input files were made in "{directory_name}".')

src/phelel/velph/cli/el_bands/plot.py

@@ -0,0 +1,113 @@
+"""Implementation of velph-el_bands-plot."""
+
+from __future__ import annotations
+
+import pathlib
+from typing import Optional
+
+import click
+import h5py
+import matplotlib.pyplot as plt
+
+from phelel.velph.cli.utils import (
+    get_distances_along_BZ_path,
+    get_reclat_from_vaspout,
+    get_special_points,
+)
+
+
+def plot_el_bandstructures(
+    window: tuple[float, float],
+    vaspout_filename_bands: pathlib.Path,
+    vaspout_filename_dos: Optional[pathlib.Path] = None,
+    plot_filename: pathlib.Path = pathlib.Path("el_bands/el_bands.pdf"),
+):
+    """Plot electronic band structure.
+
+    Parameters
+    ----------
+    window : tuple(float, float)
+        Energy window to draw band structure in eV with respect to Fermi level.
+    vaspout_filename_bands : pathlib.Path
+        Filename of vaspout.h5 of band structure.
+    vaspout_filename_dos : pathlib.Path
+        Filename of vaspout.h5 of DOS.
+    plot_filename : pathlib.Path, optional
+        File name of band structure plot.
+
+    """
+    f_h5py_bands = h5py.File(vaspout_filename_bands)
+
+    if vaspout_filename_dos:
+        f_h5py_dos = h5py.File(vaspout_filename_dos)
+        _, axs = plt.subplots(1, 2, gridspec_kw={"width_ratios": [3, 1]})
+        _plot_bands(axs[0], f_h5py_bands, window)
+        _plot_dos(axs[1], f_h5py_dos, axs[0])
+    else:
+        _, ax = plt.subplots()
+        _plot_bands(ax, f_h5py_bands, window)
+    plt.rcParams["pdf.fonttype"] = 42
+    plt.tight_layout()
+    plt.savefig(plot_filename)
+    plt.close()
+
+    click.echo(f'Electronic band structure plot was saved in "{plot_filename}".')
+
+
+def _plot_bands(ax, f_h5py, window: tuple[float, float]):
+    emin = window[0]
+    emax = window[1]
+
+    efermi = f_h5py["results"]["electron_dos"]["efermi"][()]
+    eigvals = f_h5py["results"]["electron_eigenvalues_kpoints_opt"]["eigenvalues"]
+
+    reclat = get_reclat_from_vaspout(f_h5py)
+    # k-points in reduced coordinates
+    kpoint_coords = f_h5py["results"]["electron_eigenvalues_kpoints_opt"][
+        "kpoint_coords"
+    ]
+    # Special point labels
+    labels = [
+        label.decode("utf-8")
+        for label in f_h5py["input"]["kpoints_opt"]["labels_kpoints"][:]
+    ]
+    nk_per_seg = f_h5py["input"]["kpoints_opt"]["number_kpoints"][()]
+    nk_total = len(kpoint_coords)
+    k_cart = kpoint_coords @ reclat
+    n_segments = nk_total // nk_per_seg
+    assert n_segments * nk_per_seg == nk_total
+    distances = get_distances_along_BZ_path(nk_total, n_segments, nk_per_seg, k_cart)
+    points, labels_at_points = get_special_points(
+        labels, distances, n_segments, nk_per_seg, nk_total
+    )
+
+    ax.plot(distances, eigvals[0, :, :], ".k")
+    ax.hlines(efermi, distances[0], distances[-1], "r")
+    for x in points[1:-1]:
+        ax.vlines(x, efermi + emin, efermi + emax, "k")
+    ax.set_xlim(distances[0], distances[-1])
+    ax.set_xticks(points)
+    ax.set_xticklabels(labels_at_points)
+    ax.set_ylim(efermi + emin, efermi + emax)
+    ax.set_ylabel("E[eV]", fontsize=14)
+
+
+def _plot_dos(ax, f_h5py, ax_bands):
+    efermi = f_h5py["results"]["electron_dos_kpoints_opt"]["efermi"][()]
+    dos = f_h5py["results"]["electron_dos_kpoints_opt"]["dos"][0, :]
+    energies = f_h5py["results"]["electron_dos_kpoints_opt"]["energies"][:]
+    ax.plot(dos, energies, "-k")
+    ymin, ymax = ax_bands.get_ylim()
+    i_min = 0
+    i_max = 0
+    for i, val in enumerate(energies):
+        if val < ymin:
+            i_min = i
+        if val > ymax:
+            i_max = i
+            break
+    xmax = dos[i_min:i_max].max() * 1.1
+    ax.hlines(efermi, 0, xmax, "r")
+    ax.set_xlim(0, xmax)
+    ax.set_ylim(ymin, ymax)
+    ax.yaxis.tick_right()