diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 7bd1141..f2b7a43 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -10,7 +10,7 @@ repos: - id: check-added-large-files - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.8.1 + rev: v0.8.4 hooks: - id: ruff args: [ "--fix", "--show-fixes" ] diff --git a/doc/changelog.md b/doc/changelog.md index 45f0ff9..7432987 100644 --- a/doc/changelog.md +++ b/doc/changelog.md @@ -2,6 +2,12 @@ # Change Log +## Dec-31-2024: Version 0.7.0 + +- Maintenance release to follow the change of phonopy. +- `elph_selfen_band_stop` is estimated from el-DOS instead of `elph_nbands` for + transport mode in velph. + ## Dec-9-2024: Version 0.6.6 - Collection of minor updates of velph command. diff --git a/doc/conf.py b/doc/conf.py index 1f7812a..b30883e 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -10,8 +10,8 @@ copyright = "2024, Atsushi Togo" author = "Atsushi Togo" -version = "0.6" -release = "0.6.5" +version = "0.7" +release = "0.7.0" # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration 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/pyproject.toml b/pyproject.toml index d6d01f2..a9bb5cb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,6 +12,7 @@ authors = [ ] requires-python = ">=3.9" dependencies = [ + "phonopy >= 2.33.4", "phono3py >= 3.6.0", "finufft", "click", diff --git a/src/phelel/api_phelel.py b/src/phelel/api_phelel.py index 2a4df81..35b935d 100644 --- a/src/phelel/api_phelel.py +++ b/src/phelel/api_phelel.py @@ -495,20 +495,13 @@ def _get_phonopy( supercell_matrix: Optional[Union[int, float, Sequence, np.ndarray]], primitive_matrix: Optional[Union[str, Sequence, np.ndarray]], ) -> Phonopy: - """Return Phonopy instance. - - Note - ---- - store_dense_svecs=False is necessary to be compatible with code in VASP. - - """ + """Return Phonopy instance.""" return Phonopy( self._unitcell, supercell_matrix=supercell_matrix, primitive_matrix=primitive_matrix, symprec=self._symprec, is_symmetry=self._is_symmetry, - store_dense_svecs=False, calculator=self._calculator, log_level=self._log_level, ) diff --git a/src/phelel/base/__init__.py b/src/phelel/base/__init__.py deleted file mode 100644 index 98bb25f..0000000 --- a/src/phelel/base/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Routines to calculate derivatives.""" diff --git a/src/phelel/cui/__init__.py b/src/phelel/cui/__init__.py deleted file mode 100644 index 2e4ab3c..0000000 --- a/src/phelel/cui/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Routines used for CUI.""" diff --git a/src/phelel/file_IO.py b/src/phelel/file_IO.py index bdf37f6..3dbe457 100644 --- a/src/phelel/file_IO.py +++ b/src/phelel/file_IO.py @@ -10,7 +10,7 @@ import h5py import numpy as np from phonopy.structure.atoms import PhonopyAtoms, atom_data -from phonopy.structure.cells import Primitive +from phonopy.structure.cells import Primitive, dense_to_sparse_svecs from phonopy.structure.symmetry import Symmetry from phelel.base.Dij_qij import DDijQij @@ -244,7 +244,7 @@ def _add_datasets( w.create_dataset( "primitive_masses", data=np.array(primitive.masses, dtype="double") ) - p2s_vectors, p2s_multiplicities = primitive.get_smallest_vectors() + p2s_vectors, p2s_multiplicities = _get_smallest_vectors(primitive) w.create_dataset("p2s_map", data=np.array(primitive.p2s_map, dtype="int_")) w.create_dataset("s2p_map", data=np.array(primitive.s2p_map, dtype="int_")) w.create_dataset("shortest_vectors", data=np.array(p2s_vectors, dtype="double")) @@ -308,7 +308,7 @@ def _add_datasets( data=np.array(phonon_supercell_matrix, dtype="int_", order="C"), ) if phonon_primitive is not None: - p2s_vectors, p2s_multiplicities = phonon_primitive.get_smallest_vectors() + p2s_vectors, p2s_multiplicities = _get_smallest_vectors(phonon_primitive) w.create_dataset( "phonon_p2s_map", data=np.array(phonon_primitive.p2s_map, dtype="int_") ) @@ -375,3 +375,11 @@ def _add_datasets( data=int(symmetry_dataset.uni_number), dtype="int_", ) + + +def _get_smallest_vectors(primitive: Primitive) -> tuple[np.ndarray, np.ndarray]: + """Get smallest vectors.""" + svecs, multi = primitive.get_smallest_vectors() + if primitive.store_dense_svecs: + svecs, multi = dense_to_sparse_svecs(svecs, multi) + return svecs, multi diff --git a/src/phelel/interface/__init__.py b/src/phelel/interface/__init__.py deleted file mode 100644 index 2e43ee1..0000000 --- a/src/phelel/interface/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Calculator interfaces.""" diff --git a/src/phelel/interface/vasp/__init__.py b/src/phelel/interface/vasp/__init__.py deleted file mode 100644 index 6a29c42..0000000 --- a/src/phelel/interface/vasp/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""VASP interface routines.""" diff --git a/src/phelel/utils/__init__.py b/src/phelel/utils/__init__.py deleted file mode 100644 index 9a710ae..0000000 --- a/src/phelel/utils/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Tools to calculate properties conveniently.""" 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..303320f 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]"] @@ -730,7 +757,7 @@ def _get_toml_lines( lines.append(f'spacegroup_type = "{spg_type}"') lines.append(f"tolerance = {vip.tolerance}") if len(unitcell) != len(primitive): - pmat = (np.linalg.inv(unitcell.cell) @ primitive.cell).T + pmat = (primitive.cell @ np.linalg.inv(unitcell.cell)).T lines.append("primitive_matrix = [") for v in pmat: lines.append(f" [ {v[0]:18.15f}, {v[1]:18.15f}, {v[2]:18.15f} ],") @@ -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/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/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/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..6fc89aa --- /dev/null +++ b/src/phelel/velph/cli/phonopy/generate.py @@ -0,0 +1,30 @@ +"""Implementation of velph-phonopy-generate.""" + +from __future__ import annotations + +import pathlib + +import click +import phonopy +import tomli + +from phelel.velph.cli.phelel.generate import ( + write_supercells, +) + + +def write_supercell_input_files( + toml_filename: pathlib.Path, + phonopy_yaml_filename: pathlib.Path, +) -> None: + """Generate supercells.""" + 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 + + ph = phonopy.load(phonopy_yaml_filename) + with open(toml_filename, "rb") as f: + toml_dict = tomli.load(f) + + write_supercells(ph, toml_dict, dir_name="phonopy") 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)}") diff --git a/src/phelel/velph/cli/selfenergy/generate.py b/src/phelel/velph/cli/selfenergy/generate.py index e0c508a..ee7f53e 100644 --- a/src/phelel/velph/cli/selfenergy/generate.py +++ b/src/phelel/velph/cli/selfenergy/generate.py @@ -37,7 +37,7 @@ def write_selfenergy_input_files( current_directory : Path Used for test. energy_threshold : float - Energy threshold (gap) to estimate elph_nbands used for + Energy threshold (gap) to estimate elph_selfen_band_stop used for calc_type="transport". """ @@ -83,13 +83,13 @@ def write_selfenergy_input_files( toml_incar_dict["nelm"] = 0 toml_incar_dict["elph_run"] = False - # Automatic elph_nbands setting for transport. + # Automatic elph_selfen_band_stop setting for transport. if calc_type == "transport": - # Here toml_incar_dict is updated for setting elph_nbands - band_index = _find_elph_nbands(current_directory, energy_threshold) + # Here toml_incar_dict is updated for setting elph_selfen_band_stop + band_index = _find_elph_selfen_band_stop(current_directory, energy_threshold) if band_index is not None: - click.echo(f' "elph_nbands={band_index + 1}" in INCAR is set.') - toml_incar_dict["elph_nbands"] = band_index + 1 + click.echo(f' "elph_selfen_band_stop={band_index + 1}" in INCAR is set.') + toml_incar_dict["elph_selfen_band_stop"] = band_index + 1 # POSCAR primitive = parse_cell_dict(toml_dict["primitive_cell"]) @@ -137,17 +137,19 @@ def write_selfenergy_input_files( click.echo(f'VASP input files were generated in "{directory_path}".') -def _find_elph_nbands(current_directory: pathlib.Path, energy_threshold: float) -> int: +def _find_elph_selfen_band_stop( + current_directory: pathlib.Path, energy_threshold: float +) -> int: dos_directory = current_directory / "el_bands" / "dos" if dos_directory.exists(): - click.echo('Found "el_bands/dos" directory. Estimate elph_nbands.') + click.echo('Found "el_bands/dos" directory. Estimate elph_selfen_band_stop.') vaspout_path = dos_directory / "vaspout.h5" if vaspout_path.exists(): - possiblly_occupied_band_index = _estimate_elph_nbands( + possiblly_occupied_band_index = _estimate_elph_selfen_band_stop( vaspout_path, energy_threshold=energy_threshold ) if possiblly_occupied_band_index is None: - click.echo("Estimation of elph_nbands failed.") + click.echo("Estimation of elph_selfen_band_stop failed.") return None return possiblly_occupied_band_index else: @@ -155,26 +157,23 @@ def _find_elph_nbands(current_directory: pathlib.Path, energy_threshold: float) return None -def _estimate_elph_nbands( +def _estimate_elph_selfen_band_stop( vaspout_path: pathlib.Path, energy_threshold: float = 0.5, - occupation_condition: float = 1e-5, - unoccupation_condition: float = 1e-10, + occupation_condition: float = 1e-10, ) -> Optional[int]: - """Estimate elph_nbands from eigenvalues in electronic DOS calculation. + """Estimate elph_selfen_band_stop from eigenvalues in el-DOS result. Parameters ---------- vaspout_path : pathlib.Path "vaspout.h5" path. energy_threshold : float - Energy threshold (gap) in eV to estimate elph_nbands. Default is 0.5 eV. + Energy threshold (gap) in eV to estimate elph_selfen_band_stop. Default + is 0.5 eV. occupation_condition : float Condition to determine either if bands are occupied or not. This value is used as occupation_number > occupation_condition. Default is 1e-5. - unoccupation_condition : float - Condition to determine either if bands are unoccupied or not. This value - is used as occupation_number < unoccupation_condition. Default is 1e-10. Returns ------- @@ -193,7 +192,7 @@ def _estimate_elph_nbands( nbands = eigenvalues.shape[2] unoccupied_band_index = None for i in range(nbands): - if (occupations[:, :, i] < unoccupation_condition).all(): + if (occupations[:, :, i] < occupation_condition).all(): unoccupied_band_index = i break if unoccupied_band_index is None: @@ -201,10 +200,11 @@ def _estimate_elph_nbands( occupied_eigvals = np.sort( np.extract( - occupations[:, :, unoccupied_band_index - 1] > occupation_condition, + occupations[:, :, unoccupied_band_index - 1] > occupation_condition * 0.9, eigenvalues[:, :, unoccupied_band_index - 1], ) ) + max_occupied_eigval = np.max(occupied_eigvals) for band_index in range(unoccupied_band_index, nbands): 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/src/phelel/velph/templates/__init__.py b/src/phelel/velph/templates/__init__.py index a3bf95a..36ef985 100644 --- a/src/phelel/velph/templates/__init__.py +++ b/src/phelel/velph/templates/__init__.py @@ -73,9 +73,6 @@ "elph_selfen_dw": True, "elph_selfen_delta": 0.01, "elph_selfen_temps": [300], - "elph_wf_cache_mb": 1000, - "elph_wf_cache_prefill": True, - "elph_wf_redistribute": True, "elph_nbands": 100, "elph_nbands_sum": [50, 100], "elph_selfen_gaps": True, @@ -85,7 +82,7 @@ }, "transport": { "incar": { - "elph_fermi_nedos": 5001, + "elph_fermi_nedos": 501, "elph_ismear": -24, "elph_mode": "transport", "elph_selfen_carrier_den": 0.0, @@ -144,9 +141,6 @@ 1400, ], "elph_transport_nedos": 501, - "elph_wf_redistribute": True, - "elph_wf_redistribute_opt": 0, - "elph_wf_comm_opt": 0, }, }, "el_bands": { diff --git a/src/phelel/version.py b/src/phelel/version.py index 1302b66..18ec2e8 100644 --- a/src/phelel/version.py +++ b/src/phelel/version.py @@ -1,3 +1,3 @@ """Version number.""" -__version__ = "0.6.5" +__version__ = "0.7.0" diff --git a/test/test_api_phelel.py b/test/test_api_phelel.py index e104876..3b76212 100644 --- a/test/test_api_phelel.py +++ b/test/test_api_phelel.py @@ -6,7 +6,7 @@ import numpy as np from phelel import Phelel -from phelel.file_IO import read_phelel_params_hdf5 +from phelel.file_IO import _get_smallest_vectors, read_phelel_params_hdf5 from phelel.utils.data import cmplx2real cwd = pathlib.Path(__file__).parent @@ -62,7 +62,7 @@ def _compare(filename: pathlib.Path, phe: Phelel): shortest_vectors_ref = f["shortest_vectors"][:] multiplicities_ref = f["shortest_vector_multiplicities"][:] - shortest_vectors, multiplicities = phe.primitive.get_smallest_vectors() + shortest_vectors, multiplicities = _get_smallest_vectors(phe.primitive) np.testing.assert_array_equal( shortest_vectors.shape, shortest_vectors_ref.shape ) 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/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) 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)