scnc
/
phono3py
mirror of https://github.com/phonopy/phono3py.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Refactoring

This commit is contained in:
Atsushi Togo 2025-07-19 17:11:28 +09:00
parent 35dd083b3a
commit 096269a8b6
12 changed files with 337 additions and 316 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
doc/changelog.md
View File

@ -2,6 +2,10 @@
# Change Log
## Jul-xx-2025: Version 3.17.2
- Changed `Phono3py.run_imag_self_energy()` to return `ImagSelfEnergyValues`.
## Jul-5-2025: Version 3.17.1
- Fix direct-solution crashing when executed via command line

8
doc/command-options.md
View File

@ -758,8 +758,8 @@ where the averaged phonon-phonon interaction that is read from
`kappa-mxxx(-sx-sdx).hdf5` file is used if it exists in the file. Therefore the
averaged phonon-phonon interaction has to be stored before using this option
(see {ref}`--full-pp <full_pp_option>`). The calculation result **overwrites**
`kappa-mxxx(-sx-sdx).hdf5` file. Therefore to use this option together with `-o`
option is strongly recommended.
`kappa-mxxx(-sx-sdx).hdf5` file. Therefore the original
`kappa-mxxx(-sx-sdx).hdf5` file should be backed up.
First, run full conductivity calculation,
@ -1224,10 +1224,6 @@ set. Other filters (`lzf` or integer values of 0 to 9) may be used, see h5py
documentation
(<http://docs.h5py.org/en/stable/high/dataset.html#filter-pipeline>).
### `-o`, `-i`, `--io`
These options are deprecated.
<!-- (output_filename_option)=
### `-o` (command option only)

125
phono3py/api_phono3py.py
View File

@ -36,6 +36,8 @@
from __future__ import annotations
import copy
import dataclasses
import os
import warnings
from collections.abc import Sequence
from typing import Literal, cast
@ -107,6 +109,16 @@ from phono3py.phonon3.spectral_function import run_spectral_function
from phono3py.version import __version__
@dataclasses.dataclass
class ImagSelfEnergyValues:
"""Parameters for imaginary self-energy calculation."""
frequency_points: NDArray | None
gammas: NDArray
scattering_event_class: int | None = None
detailed_gammas: Sequence | None = None
class Phono3py:
"""Phono3py main class.
@ -152,7 +164,7 @@ class Phono3py:
self,
unitcell: PhonopyAtoms,
supercell_matrix: ArrayLike | None = None,
primitive_matrix: ArrayLike | None = None,
primitive_matrix: str | ArrayLike | None = None,
phonon_supercell_matrix: ArrayLike | None = None,
cutoff_frequency: float = 1e-4,
frequency_factor_to_THz: float | None = None,
@ -244,8 +256,7 @@ class Phono3py:
self._supercell_matrix = np.array(
shape_supercell_matrix(supercell_matrix), dtype="int64", order="C"
)
pmat = self._determine_primitive_matrix(primitive_matrix)
self._primitive_matrix = pmat
self._primitive_matrix = self._determine_primitive_matrix(primitive_matrix)
self._nac_params = None
if phonon_supercell_matrix is not None:
self._phonon_supercell_matrix = np.array(
@ -279,18 +290,17 @@ class Phono3py:
self._search_phonon_supercell_symmetry()
# Displacements and supercells
self._supercells_with_displacements = None
self._dataset = None
self._phonon_dataset = None
self._phonon_supercells_with_displacements = None
self._dataset: dict | None = None
self._phonon_dataset: dict | None = None
self._supercells_with_displacements: list | None = None
self._phonon_supercells_with_displacements: list | None = None
# Thermal conductivity
# conductivity_RTA or conductivity_LBTE class instance
self._thermal_conductivity = None
# Imaginary part of self energy at frequency points
self._gammas = None
self._scattering_event_class = None
self._ise_params = None
# Frequency shift (real part of bubble diagram)
self._real_self_energy = None
@ -307,9 +317,9 @@ class Phono3py:
# MLP
self._mlp = None
self._mlp_dataset = None
self._mlp_dataset: dict | None = None
self._phonon_mlp = None
self._phonon_mlp_dataset = None
self._phonon_mlp_dataset: dict | None = None
# Setup interaction
self._interaction = None
@ -584,10 +594,10 @@ class Phono3py:
return self._phonon_supercell_matrix
@property
def primitive_matrix(self) -> NDArray:
def primitive_matrix(self) -> NDArray | None:
"""Return transformation matrix to primitive cell from unit cell.
ndarray
ndarray or None
Primitive matrix with respect to unit cell.
shape=(3, 3), dtype='double', order='C'
@ -805,7 +815,7 @@ class Phono3py:
self._phonon_supercell.masses = s_masses
@property
def supercells_with_displacements(self) -> list[PhonopyAtoms]:
def supercells_with_displacements(self) -> list[PhonopyAtoms | None]:
"""Return supercells with displacements.
list of PhonopyAtoms
@ -815,6 +825,7 @@ class Phono3py:
"""
if self._supercells_with_displacements is None:
self._build_supercells_with_displacements()
assert self._supercells_with_displacements is not None
return self._supercells_with_displacements
@property
@ -884,7 +895,7 @@ class Phono3py:
"""
dataset = self._dataset
if self._dataset is None:
if dataset is None:
raise RuntimeError("displacement dataset is not set.")
if "first_atoms" in dataset:
@ -1079,7 +1090,9 @@ class Phono3py:
@property
def detailed_gammas(self):
"""Return detailed gamma."""
return self._detailed_gammas
if self._ise_params is None:
raise RuntimeError("Imaginary self energy parameters are not set.")
return self._ise_params.detailed_gammas
@property
def grid(self):
@ -1093,12 +1106,12 @@ class Phono3py:
def init_phph_interaction(
self,
nac_q_direction=None,
constant_averaged_interaction=None,
frequency_scale_factor=None,
nac_q_direction: ArrayLike | None = None,
constant_averaged_interaction: float | None = None,
frequency_scale_factor: float | None = None,
symmetrize_fc3q: bool = False,
lapack_zheev_uplo="L",
openmp_per_triplets=None,
lapack_zheev_uplo: Literal["L", "U"] = "L",
openmp_per_triplets: bool | None = None,
):
"""Initialize ph-ph interaction calculation.
@ -1252,7 +1265,7 @@ class Phono3py:
self,
distance: float | None = None,
cutoff_pair_distance: float | None = None,
is_plusminus: bool | str = "auto",
is_plusminus: bool | Literal["auto"] = "auto",
is_diagonal: bool = True,
number_of_snapshots: int | Literal["auto"] | None = None,
random_seed: int | None = None,
@ -1398,7 +1411,7 @@ class Phono3py:
def generate_fc2_displacements(
self,
distance: float | None = None,
is_plusminus: str = "auto",
is_plusminus: bool | Literal["auto"] = "auto",
is_diagonal: bool = False,
number_of_snapshots: int | Literal["auto"] | None = None,
random_seed: int | None = None,
@ -1811,7 +1824,7 @@ class Phono3py:
write_gamma_detail=False,
keep_gamma_detail=False,
output_filename=None,
):
) -> ImagSelfEnergyValues:
"""Calculate imaginary part of self-energy of bubble diagram (Gamma).
Pi = Delta - i Gamma.
@ -1876,7 +1889,6 @@ class Phono3py:
else:
self._temperatures = temperatures
self._grid_points = grid_points
self._scattering_event_class = scattering_event_class
vals = get_imag_self_energy(
self._interaction,
grid_points,
@ -1894,25 +1906,34 @@ class Phono3py:
log_level=self._log_level,
)
if keep_gamma_detail:
(self._frequency_points, self._gammas, self._detailed_gammas) = vals
self._ise_params = ImagSelfEnergyValues(
frequency_points=vals[0],
gammas=vals[1],
scattering_event_class=scattering_event_class,
detailed_gammas=vals[2],
)
else:
self._frequency_points, self._gammas = vals
self._ise_params = ImagSelfEnergyValues(
frequency_points=vals[0], gammas=vals[1]
)
if write_txt:
self._write_imag_self_energy(output_filename=output_filename)
return vals
return self._ise_params
def _write_imag_self_energy(self, output_filename=None):
if self._ise_params is None:
raise RuntimeError("Imaginary self-energy is not calculated.")
write_imag_self_energy(
self._gammas,
self._ise_params.gammas,
self.mesh_numbers,
self._grid_points,
self._band_indices,
self._frequency_points,
self._ise_params.frequency_points,
self._temperatures,
self._sigmas,
scattering_event_class=self._scattering_event_class,
scattering_event_class=self._ise_params.scattering_event_class,
output_filename=output_filename,
is_mesh_symmetry=self._is_mesh_symmetry,
log_level=self._log_level,
@ -2352,7 +2373,9 @@ class Phono3py:
)
def save(
self, filename: str = "phono3py_params.yaml", settings: dict | None = None
self,
filename: str | os.PathLike = "phono3py_params.yaml",
settings: dict | None = None,
):
"""Save parameters in Phono3py instants into file.
@ -2470,9 +2493,9 @@ class Phono3py:
test_size=test_size,
)
def save_phonon_mlp(self, filename: str | None = None):
def save_phonon_mlp(self, filename: str | os.PathLike | None = None):
"""Save machine learning potential."""
if self._mlp is None:
if self._phonon_mlp is None:
raise RuntimeError("MLP is not developed yet.")
self._phonon_mlp.save(filename=filename)
@ -2495,8 +2518,10 @@ class Phono3py:
constants.
"""
if self._mlp is None and self._phonon_mlp is None:
if self._mlp is None:
raise RuntimeError("MLP is not developed yet.")
if self._phonon_mlp is None:
raise RuntimeError("Phonon MLP is not developed yet.")
if self.phonon_supercells_with_displacements is None:
raise RuntimeError(
@ -2633,15 +2658,22 @@ class Phono3py:
return supercells
def _build_supercells_with_displacements(self):
assert self._dataset is not None
magmoms = self._supercell.magnetic_moments
masses = self._supercell.masses
numbers = self._supercell.numbers
lattice = self._supercell.cell
supercells = self._build_phonon_supercells_with_displacements(
self._supercell, self._dataset
# One displacement supercells
supercells = cast(
list[PhonopyAtoms | None],
self._build_phonon_supercells_with_displacements(
self._supercell, self._dataset
),
)
# Two displacement supercells
if "first_atoms" in self._dataset:
for disp1 in self._dataset["first_atoms"]:
disp_cart1 = disp1["displacement"]
@ -2679,7 +2711,9 @@ class Phono3py:
return get_primitive(supercell, t_mat, self._symprec, store_dense_svecs=True)
def _determine_primitive_matrix(self, primitive_matrix):
def _determine_primitive_matrix(
self, primitive_matrix: str | ArrayLike | None
) -> NDArray | None:
pmat = get_primitive_matrix(primitive_matrix, symprec=self._symprec)
if isinstance(pmat, str) and pmat == "auto":
return guess_primitive_matrix(self._unitcell, symprec=self._symprec)
@ -2719,7 +2753,8 @@ class Phono3py:
nac_params=self._nac_params,
)
freqs, _, _ = self.get_phonon_data()
gp_Gamma = self._bz_grid.gp_Gamma
gp_Gamma = self._interaction.bz_grid.gp_Gamma
assert freqs is not None
if np.sum(freqs[gp_Gamma] < self._cutoff_frequency) < 3:
for i, f in enumerate(freqs[gp_Gamma, :3]):
if not (f < self._cutoff_frequency):
@ -2775,6 +2810,9 @@ class Phono3py:
def _set_forces_energies(
self, values, target: Literal["forces", "supercell_energies"]
):
if self._dataset is None:
raise RuntimeError("Dataset is not available.")
if "first_atoms" in self._dataset: # type-1
count = 0
for disp1 in self._dataset["first_atoms"]:
@ -2858,13 +2896,11 @@ class Phono3py:
is_plusminus: bool = False,
random_seed: int | None = None,
max_distance: float | None = None,
):
) -> dict:
if random_seed is not None and random_seed >= 0 and random_seed < 2**32:
_random_seed = random_seed
dataset = {"random_seed": _random_seed}
else:
_random_seed = None
dataset = {}
d = get_random_displacements_dataset(
number_of_snapshots,
number_of_atoms,
@ -2873,7 +2909,10 @@ class Phono3py:
is_plusminus=is_plusminus,
max_distance=max_distance,
)
dataset["displacements"] = d
if _random_seed is None:
dataset = {"displacements": d}
else:
dataset = {"random_seed": _random_seed, "displacements": d}
return dataset
def _check_mlp_dataset(self, mlp_dataset: dict):

116
phono3py/cui/create_force_constants.py
View File

@ -183,18 +183,9 @@ def _read_FORCES_FC3_or_FC2(
return dataset
def run_pypolymlp_to_compute_forces(
def develop_pypolymlp(
ph3py: Phono3py,
mlp_params: str | dict | PypolymlpParams | None = None,
displacement_distance: float | None = None,
number_of_snapshots: int | Literal["auto"] | None = None,
number_estimation_factor: int | None = None,
random_seed: int | None = None,
prepare_dataset: bool = False,
fc_calculator: str | None = None,
fc_calculator_options: str | None = None,
cutoff_pair_distance: float | None = None,
symfc_memory_size: float | None = None,
mlp_filename: str | os.PathLike | None = None,
log_level: int = 0,
):
@ -253,58 +244,71 @@ def run_pypolymlp_to_compute_forces(
if log_level:
print("-" * 30 + " pypolymlp end " + "-" * 31, flush=True)
if prepare_dataset:
if displacement_distance is None:
_displacement_distance = 0.01
def generate_displacements_and_evaluate_pypolymlp(
ph3py: Phono3py,
displacement_distance: float | None = None,
number_of_snapshots: int | Literal["auto"] | None = None,
number_estimation_factor: int | None = None,
random_seed: int | None = None,
fc_calculator: str | None = None,
fc_calculator_options: str | None = None,
cutoff_pair_distance: float | None = None,
symfc_memory_size: float | None = None,
log_level: int = 0,
):
"""Generate displacements and evaluate forces by pypolymlp."""
if displacement_distance is None:
_displacement_distance = 0.01
else:
_displacement_distance = displacement_distance
if log_level:
if number_of_snapshots:
print("Generate random displacements")
print(
" Twice of number of snapshots will be generated "
"for plus-minus displacements."
)
else:
_displacement_distance = displacement_distance
if log_level:
if number_of_snapshots:
print("Generate random displacements")
print(
" Twice of number of snapshots will be generated "
"for plus-minus displacements."
)
else:
print("Generate displacements")
print(
f" Displacement distance: {_displacement_distance:.5f}".rstrip(
"0"
).rstrip(".")
print("Generate displacements")
print(
f" Displacement distance: {_displacement_distance:.5f}".rstrip("0").rstrip(
"."
)
cutoff_pair_distance = determine_cutoff_pair_distance(
fc_calculator=fc_calculator,
fc_calculator_options=fc_calculator_options,
cutoff_pair_distance=cutoff_pair_distance,
symfc_memory_size=symfc_memory_size,
random_displacements=number_of_snapshots,
supercell=ph3py.supercell,
primitive=ph3py.primitive,
symmetry=ph3py.symmetry,
log_level=log_level,
)
ph3py.generate_displacements(
distance=_displacement_distance,
cutoff_pair_distance=cutoff_pair_distance,
is_plusminus=True,
number_of_snapshots=number_of_snapshots,
random_seed=random_seed,
number_estimation_factor=number_estimation_factor,
)
if log_level:
print(
f"Evaluate forces in {ph3py.displacements.shape[0]} supercells "
"by pypolymlp",
flush=True,
)
cutoff_pair_distance = determine_cutoff_pair_distance(
fc_calculator=fc_calculator,
fc_calculator_options=fc_calculator_options,
cutoff_pair_distance=cutoff_pair_distance,
symfc_memory_size=symfc_memory_size,
random_displacements=number_of_snapshots,
supercell=ph3py.supercell,
primitive=ph3py.primitive,
symmetry=ph3py.symmetry,
log_level=log_level,
)
ph3py.generate_displacements(
distance=_displacement_distance,
cutoff_pair_distance=cutoff_pair_distance,
is_plusminus=True,
number_of_snapshots=number_of_snapshots,
random_seed=random_seed,
number_estimation_factor=number_estimation_factor,
)
if ph3py.supercells_with_displacements is None:
raise RuntimeError("Displacements are not set. Run generate_displacements.")
if log_level:
print(
f"Evaluate forces in {ph3py.displacements.shape[0]} supercells "
"by pypolymlp",
flush=True,
)
ph3py.evaluate_mlp()
if ph3py.supercells_with_displacements is None:
raise RuntimeError("Displacements are not set. Run generate_displacements.")
ph3py.evaluate_mlp()
def run_pypolymlp_to_compute_phonon_forces(

6
phono3py/cui/create_force_sets.py
View File

@ -200,17 +200,15 @@ def create_FORCES_FC2_from_FORCE_SETS(log_level):
def create_FORCE_SETS_from_FORCES_FCx(
phonon_smat, input_filename: Optional[str], cell_filename: Optional[str], log_level
phonon_smat, cell_filename: Optional[str], log_level
):
"""Convert FORCES_FC3 or FORCES_FC2 to FORCE_SETS."""
if cell_filename is not None and is_file_phonopy_yaml(
cell_filename, keyword="phono3py"
):
disp_filename = cell_filename
elif input_filename is None:
disp_filename = "phono3py_disp.yaml"
else:
disp_filename = f"phono3py_disp.{input_filename}.yaml"
disp_filename = "phono3py_disp.yaml"
if phonon_smat is not None:
forces_filename = "FORCES_FC2"
else:

6
phono3py/cui/load.py
View File

@ -51,8 +51,8 @@ from phonopy.structure.cells import determinant
from phono3py import Phono3py
from phono3py.cui.create_force_constants import (
develop_pypolymlp,
parse_forces,
run_pypolymlp_to_compute_forces,
)
from phono3py.file_IO import read_fc2_from_hdf5, read_fc3_from_hdf5
from phono3py.interface.fc_calculator import (
@ -371,9 +371,7 @@ def load(
if produce_fc:
if ph3py.fc3 is None and use_pypolymlp:
run_pypolymlp_to_compute_forces(
ph3py, mlp_params=mlp_params, log_level=log_level
)
develop_pypolymlp(ph3py, mlp_params=mlp_params, log_level=log_level)
compute_force_constants_from_datasets(
ph3py,

46
phono3py/cui/phono3py_argparse.py
View File

@ -58,7 +58,7 @@ def get_parser(load_phono3py_yaml: bool = False):
"--alm",
dest="use_alm",
action="store_true",
default=False,
default=None,
help=("Use ALM for generating 2nd and 3rd force constants in one fitting"),
)
parser.add_argument(
@ -180,7 +180,7 @@ def get_parser(load_phono3py_yaml: bool = False):
"--create-force-sets",
dest="force_sets_mode",
action="store_true",
default=False,
default=None,
help="Create phonopy FORCE_SETS from FORCES_FC2",
)
parser.add_argument(
@ -188,7 +188,7 @@ def get_parser(load_phono3py_yaml: bool = False):
"--collective-phonon",
dest="solve_collective_phonon",
action="store_true",
default=False,
default=None,
help="Solve collective phonons",
)
if load_phono3py_yaml:
@ -280,14 +280,14 @@ def get_parser(load_phono3py_yaml: bool = False):
"--fc2",
dest="read_fc2",
action="store_true",
default=False,
default=None,
help="Read second order force constants",
)
parser.add_argument(
"--fc3",
dest="read_fc3",
action="store_true",
default=False,
default=None,
help="Read third order force constants",
)
parser.add_argument(
@ -338,7 +338,7 @@ def get_parser(load_phono3py_yaml: bool = False):
"--fs2f2",
"--force-sets-to-forces-fc2",
dest="force_sets_to_forces_fc2_mode",
default=False,
default=None,
action="store_true",
help="Create FORCES_FC2 from FORCE_SETS",
)
@ -346,7 +346,7 @@ def get_parser(load_phono3py_yaml: bool = False):
"--full-pp",
dest="is_full_pp",
action="store_true",
default=False,
default=None,
help=(
"Calculate full ph-ph interaction for RTA conductivity."
"This may be activated when full elements of ph-ph interaction "
@ -383,14 +383,14 @@ def get_parser(load_phono3py_yaml: bool = False):
"--generalized-regular-grid",
dest="use_grg",
action="store_true",
default=False,
default=None,
help="Use generalized regular grid.",
)
parser.add_argument(
"--gruneisen",
dest="is_gruneisen",
action="store_true",
default=False,
default=None,
help="Calculate phonon Gruneisen parameter",
)
parser.add_argument(
@ -406,21 +406,11 @@ def get_parser(load_phono3py_yaml: bool = False):
default=None,
help="hdf5 compression filter (default: gzip)",
)
if not load_phono3py_yaml:
parser.add_argument(
"-i", dest="input_filename", default=None, help="Input filename extension"
)
parser.add_argument(
"--io",
dest="input_output_filename",
default=None,
help="Input and output filename extension",
)
parser.add_argument(
"--ion-clamped",
dest="ion_clamped",
action="store_true",
default=False,
default=None,
help=(
"Atoms are clamped under applied strain in Gruneisen parameter calculation"
),
@ -429,35 +419,35 @@ def get_parser(load_phono3py_yaml: bool = False):
"--ise",
dest="is_imag_self_energy",
action="store_true",
default=False,
default=None,
help="Calculate imaginary part of self energy",
)
parser.add_argument(
"--isotope",
dest="is_isotope",
action="store_true",
default=False,
default=None,
help="Isotope scattering lifetime",
)
parser.add_argument(
"--jdos",
dest="is_joint_dos",
action="store_true",
default=False,
default=None,
help="Calculate joint density of states",
)
parser.add_argument(
"--kubo",
dest="is_kubo_kappa",
action="store_true",
default=False,
default=None,
help="Choose Kubo lattice thermal conductivity.",
)
parser.add_argument(
"--lbte",
dest="is_lbte",
action="store_true",
default=False,
default=None,
help="Calculate thermal conductivity LBTE with Chaput's method",
)
parser.add_argument(
@ -546,7 +536,7 @@ def get_parser(load_phono3py_yaml: bool = False):
"--nomeshsym",
dest="is_nomeshsym",
action="store_true",
default=False,
default=None,
help="No symmetrization of triplets is made.",
)
if load_phono3py_yaml:
@ -595,10 +585,6 @@ def get_parser(load_phono3py_yaml: bool = False):
default=None,
help="Output yaml filename instead of default filename of phono3py.yaml",
)
else:
parser.add_argument(
"-o", dest="output_filename", default=None, help="Output filename extension"
)
parser.add_argument(
"--pa",
"--primitive-axis",

191
phono3py/cui/phono3py_script.py
View File

@ -40,7 +40,6 @@ import argparse
import os
import pathlib
import sys
import warnings
from collections.abc import Sequence
from typing import cast
@ -69,7 +68,10 @@ from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import isclose as cells_isclose
from phono3py import Phono3py, Phono3pyIsotope, Phono3pyJointDos
from phono3py.cui.create_force_constants import run_pypolymlp_to_compute_forces
from phono3py.cui.create_force_constants import (
develop_pypolymlp,
generate_displacements_and_evaluate_pypolymlp,
)
from phono3py.cui.create_force_sets import (
create_FORCE_SETS_from_FORCES_FCx,
create_FORCES_FC2_from_FORCE_SETS,
@ -193,7 +195,6 @@ def _finalize_phono3py(
def _get_run_mode(settings: Phono3pySettings):
"""Extract run mode from settings."""
run_mode = None
if settings.is_gruneisen:
run_mode = "gruneisen"
elif settings.is_joint_dos:
@ -216,6 +217,8 @@ def _get_run_mode(settings: Phono3pySettings):
run_mode = "displacements"
elif settings.write_phonon:
run_mode = "phonon"
else:
run_mode = "force constants"
return run_mode
@ -314,33 +317,6 @@ def _read_phono3py_settings(
return settings, confs_dict, cell_filename
def _get_input_output_filenames_from_args(args: argparse.Namespace):
"""Return strings inserted to input and output filenames."""
if args.input_filename is not None:
warnings.warn(
"-i option is deprecated and will be removed soon.",
DeprecationWarning,
stacklevel=2,
)
input_filename = args.input_filename
if args.output_filename is not None:
warnings.warn(
"-o option is deprecated and will be removed soon.",
DeprecationWarning,
stacklevel=2,
)
output_filename = args.output_filename
if args.input_output_filename is not None:
warnings.warn(
"--io option is deprecated and will be removed soon.",
DeprecationWarning,
stacklevel=2,
)
input_filename = args.input_output_filename
output_filename = args.input_output_filename
return input_filename, output_filename
def _get_default_values(settings: Phono3pySettings):
"""Set default values."""
# Brillouin zone integration: Tetrahedron (default) or smearing method
@ -531,8 +507,9 @@ def _create_supercells_with_displacements(
confs_dict: dict,
unitcell_filename: str,
interface_mode: str | None,
symprec: float,
log_level: int,
symprec: float = 1e-5,
output_yaml_filename: str | None = None,
log_level: int = 0,
):
"""Create supercells and write displacements."""
if (
@ -540,7 +517,7 @@ def _create_supercells_with_displacements(
or settings.random_displacements is not None
or settings.random_displacements_fc2 is not None
):
phono3py = create_phono3py_supercells(
ph3py = create_phono3py_supercells(
cell_info,
settings,
symprec,
@ -550,7 +527,7 @@ def _create_supercells_with_displacements(
if pathlib.Path("BORN").exists():
store_nac_params(
cast(Phonopy, phono3py),
cast(Phonopy, ph3py),
cast(PhonopySettings, settings),
cell_info.phono3py_yaml,
unitcell_filename,
@ -559,20 +536,69 @@ def _create_supercells_with_displacements(
)
if log_level:
if phono3py.supercell.magnetic_moments is None:
print("Spacegroup: %s" % phono3py.symmetry.get_international_table())
if ph3py.supercell.magnetic_moments is None:
print("Spacegroup: %s" % ph3py.symmetry.get_international_table())
else:
print(
"Number of symmetry operations in supercell: %d"
% len(phono3py.symmetry.symmetry_operations["rotations"])
% len(ph3py.symmetry.symmetry_operations["rotations"])
)
ph3py.save("phono3py_disp.yaml")
if log_level > 0:
print('Displacement dataset was written in "phono3py_disp.yaml".')
_finalize_phono3py(
phono3py,
ph3py,
confs_dict,
log_level,
write_displacements=True,
filename="phono3py_disp.yaml",
filename=output_yaml_filename,
)
def _run_pypolymlp(
ph3py: Phono3py,
settings: Phono3pySettings,
confs_dict: dict,
output_yaml_filename: str | None = None,
mlp_eval_filename: str = "phono3py_mlp_eval_dataset.yaml",
log_level: int = 0,
):
assert ph3py.mlp_dataset is None
if ph3py.dataset is not None: # If None, load mlp from polymlp.yaml.
ph3py.mlp_dataset = ph3py.dataset
ph3py.dataset = None
develop_pypolymlp(
ph3py,
mlp_params=settings.mlp_params,
log_level=log_level,
)
if settings.create_displacements or settings.random_displacements is not None:
generate_displacements_and_evaluate_pypolymlp(
ph3py,
displacement_distance=settings.displacement_distance,
number_of_snapshots=settings.random_displacements,
number_estimation_factor=settings.rd_number_estimation_factor,
random_seed=settings.random_seed,
fc_calculator=settings.fc_calculator,
fc_calculator_options=settings.fc_calculator_options,
cutoff_pair_distance=settings.cutoff_pair_distance,
symfc_memory_size=settings.symfc_memory_size,
log_level=log_level,
)
if log_level:
print(f'Dataset generated using MLPs was written in "{mlp_eval_filename}".')
ph3py.save(mlp_eval_filename)
else:
if log_level:
print(
"Generate displacements (--rd or -d) for proceeding to phonon "
"calculations."
)
_finalize_phono3py(
ph3py, confs_dict, log_level=log_level, filename=output_yaml_filename
)
@ -695,7 +721,6 @@ def _produce_force_constants(
def _run_gruneisen_then_exit(
phono3py: Phono3py,
settings: Phono3pySettings,
output_filename: str | os.PathLike | None,
log_level: int,
):
"""Run mode Grueneisen parameter calculation from fc3."""
@ -742,7 +767,6 @@ def _run_gruneisen_then_exit(
ion_clamped=settings.ion_clamped,
factor=get_physical_units().DefaultToTHz,
symprec=phono3py.symmetry.tolerance,
output_filename=output_filename,
log_level=log_level,
)
@ -755,7 +779,6 @@ def _run_jdos_then_exit(
phono3py: Phono3py,
settings: Phono3pySettings,
updated_settings: dict,
output_filename: str | None,
log_level: int,
):
"""Run joint-DOS calculation."""
@ -777,7 +800,6 @@ def _run_jdos_then_exit(
use_grg=settings.use_grg,
is_mesh_symmetry=settings.is_mesh_symmetry,
symprec=phono3py.symmetry.tolerance,
output_filename=output_filename,
log_level=log_level,
)
@ -843,8 +865,6 @@ def _init_phph_interaction(
phono3py: Phono3py,
settings: Phono3pySettings,
updated_settings: dict,
input_filename: str | None,
output_filename: str | None,
log_level: int,
):
"""Initialize ph-ph interaction and phonons on grid."""
@ -901,7 +921,6 @@ def _init_phph_interaction(
ir_grid_points=ir_grid_points,
ir_grid_weights=ir_grid_weights,
compression=settings.hdf5_compression,
filename=output_filename,
)
if filename:
if log_level:
@ -913,9 +932,7 @@ def _init_phph_interaction(
sys.exit(1)
if settings.read_phonon:
phonons = read_phonon_from_hdf5(
phono3py.mesh_numbers, filename=input_filename, verbose=(log_level > 0)
)
phonons = read_phonon_from_hdf5(phono3py.mesh_numbers, verbose=(log_level > 0))
if phonons is None:
print("Reading phonons failed.")
if log_level:
@ -981,12 +998,10 @@ def main(**argparse_control):
else:
args, log_level = _start_phono3py(**argparse_control)
output_yaml_filename: str | None
if load_phono3py_yaml:
input_filename = None
output_filename = None
output_yaml_filename = args.output_yaml_filename
else:
(input_filename, output_filename) = _get_input_output_filenames_from_args(args)
output_yaml_filename = None
settings, confs_dict, cell_filename = _read_phono3py_settings(
@ -1000,7 +1015,7 @@ def main(**argparse_control):
sys.exit(0)
if args.force_sets_mode:
create_FORCE_SETS_from_FORCES_FCx(
settings.phonon_supercell_matrix, input_filename, cell_filename, log_level
settings.phonon_supercell_matrix, cell_filename, log_level
)
if log_level:
print_end_phono3py()
@ -1055,6 +1070,7 @@ def main(**argparse_control):
if run_mode is None:
run_mode = _get_run_mode(settings)
assert run_mode is not None
######################################################
# Create supercells with displacements and then exit #
@ -1066,8 +1082,9 @@ def main(**argparse_control):
confs_dict,
unitcell_filename,
interface_mode,
symprec,
log_level,
symprec=symprec,
output_yaml_filename=output_yaml_filename,
log_level=log_level,
)
#######################
@ -1091,8 +1108,6 @@ def main(**argparse_control):
run_mode,
ph3py,
unitcell_filename,
input_filename,
output_filename,
interface_mode,
)
@ -1247,49 +1262,14 @@ def main(**argparse_control):
# polynomial MLPs #
###################
if settings.use_pypolymlp:
assert ph3py.mlp_dataset is None
if ph3py.dataset is not None: # If None, load mlp from polymlp.yaml.
ph3py.mlp_dataset = ph3py.dataset
ph3py.dataset = None
prepare_dataset = (
settings.create_displacements or settings.random_displacements is not None
)
run_pypolymlp_to_compute_forces(
_run_pypolymlp(
ph3py,
mlp_params=settings.mlp_params,
displacement_distance=settings.displacement_distance,
number_of_snapshots=settings.random_displacements,
number_estimation_factor=settings.rd_number_estimation_factor,
random_seed=settings.random_seed,
fc_calculator=settings.fc_calculator,
fc_calculator_options=settings.fc_calculator_options,
cutoff_pair_distance=settings.cutoff_pair_distance,
symfc_memory_size=settings.symfc_memory_size,
prepare_dataset=prepare_dataset,
settings,
confs_dict,
output_yaml_filename=output_yaml_filename,
log_level=log_level,
)
if ph3py.dataset is not None:
mlp_eval_filename = "phono3py_mlp_eval_dataset.yaml"
if log_level:
print(
"Dataset generated using MLPs was written in "
f'"{mlp_eval_filename}".'
)
ph3py.save(mlp_eval_filename)
# pypolymlp dataset is stored in "polymlp.yaml" and stop here.
if not prepare_dataset:
if log_level:
print(
"Generate displacements (--rd or -d) for proceeding to phonon "
"calculations."
)
_finalize_phono3py(
ph3py, confs_dict, log_level, filename=output_yaml_filename
)
###########################
# Produce force constants #
###########################
@ -1299,21 +1279,19 @@ def main(**argparse_control):
# Phonon Gruneisen parameter and then exit #
############################################
if settings.is_gruneisen:
_run_gruneisen_then_exit(ph3py, settings, output_filename, log_level)
_run_gruneisen_then_exit(ph3py, settings, log_level)
#################
# Show settings #
#################
if log_level and run_mode is not None:
if log_level:
show_phono3py_settings(ph3py, settings, updated_settings, log_level)
###########################
# Joint DOS and then exit #
###########################
if run_mode == "jdos":
_run_jdos_then_exit(
ph3py, settings, updated_settings, output_filename, log_level
)
_run_jdos_then_exit(ph3py, settings, updated_settings, log_level)
################################################
# Mass variances for phonon-isotope scattering #
@ -1343,13 +1321,11 @@ def main(**argparse_control):
########################################
# Initialize phonon-phonon interaction #
########################################
if run_mode is not None:
if run_mode in run_modes_with_mesh:
_init_phph_interaction(
ph3py,
settings,
updated_settings,
input_filename,
output_filename,
log_level,
)
@ -1367,7 +1343,6 @@ def main(**argparse_control):
scattering_event_class=settings.scattering_event_class,
write_txt=True,
write_gamma_detail=settings.write_gamma_detail,
output_filename=output_filename,
)
#####################################################
@ -1382,7 +1357,6 @@ def main(**argparse_control):
num_frequency_points=updated_settings["num_frequency_points"],
write_txt=True,
write_hdf5=True,
output_filename=output_filename,
)
#######################################################
@ -1398,7 +1372,6 @@ def main(**argparse_control):
num_points_in_batch=updated_settings["num_points_in_batch"],
write_txt=True,
write_hdf5=True,
output_filename=output_filename,
)
####################################
@ -1435,8 +1408,6 @@ def main(**argparse_control):
read_pp=settings.read_pp,
write_LBTE_solution=settings.write_LBTE_solution,
compression=settings.hdf5_compression,
input_filename=input_filename,
output_filename=output_filename,
)
else:
if log_level:

21
phono3py/cui/show_log.py
View File

@ -46,13 +46,11 @@ from phono3py.cui.settings import Phono3pySettings
def show_general_settings(
settings,
run_mode,
phono3py,
cell_filename,
input_filename,
output_filename,
interface_mode,
settings: Phono3pySettings,
run_mode: str,
phono3py: Phono3py,
cell_filename: str,
interface_mode: str | None,
):
"""Show general setting information."""
is_primitive_axes_auto = (
@ -65,11 +63,10 @@ def show_general_settings(
print("-" * 29 + " General settings " + "-" * 29)
if run_mode:
print(f"Run mode: {run_mode}")
if output_filename:
print(f"Output filename is modified by {output_filename}.")
if input_filename:
print(f"Input filename is modified by {input_filename}.")
if settings.use_pypolymlp:
print(f"Run mode: pypolymlp & {run_mode}")
else:
print(f"Run mode: {run_mode}")
if settings.hdf5_compression:
print(f"HDF5 data compression filter: {settings.hdf5_compression}")
if interface_mode:

9
phono3py/phonon3/displacement_fc3.py
View File

@ -34,6 +34,10 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
from typing import Literal
import numpy as np
from phonopy.harmonic.displacement import (
directions_axis,
@ -131,7 +135,10 @@ def direction_to_displacement(
def get_third_order_displacements(
cell: PhonopyAtoms, symmetry: Symmetry, is_plusminus="auto", is_diagonal=False
cell: PhonopyAtoms,
symmetry: Symmetry,
is_plusminus: bool | Literal["auto"] = "auto",
is_diagonal: bool = False,
):
"""Create displacement dataset.

48
phono3py/phonon3/imag_self_energy.py
View File

@ -34,11 +34,14 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import sys
import warnings
from typing import List, Optional
from collections.abc import Sequence
import numpy as np
from numpy.typing import ArrayLike, NDArray
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.physical_units import get_physical_units
@ -576,20 +579,20 @@ class ImagSelfEnergy:
def get_imag_self_energy(
interaction: Interaction,
grid_points,
temperatures,
sigmas=None,
frequency_points=None,
frequency_step=None,
num_frequency_points=None,
frequency_points_at_bands=False,
num_points_in_batch=None,
scattering_event_class=None, # class 1 or 2
write_gamma_detail=False,
return_gamma_detail=False,
output_filename=None,
log_level=0,
):
grid_points: ArrayLike,
temperatures: ArrayLike,
sigmas: Sequence[float | None] | None = None,
frequency_points: ArrayLike | None = None,
frequency_step: float | None = None,
num_frequency_points: int | None = None,
frequency_points_at_bands: bool = False,
num_points_in_batch: int | None = None,
scattering_event_class: int | None = None, # class 1 or 2
write_gamma_detail: bool = False,
return_gamma_detail: bool = False,
output_filename: str | None = None,
log_level: int = 0,
) -> tuple[NDArray | None, NDArray, Sequence]:
"""Imaginary-part of self-energy at frequency points.
Band indices to be calculated at are found in Interaction instance.
@ -599,12 +602,12 @@ def get_imag_self_energy(
interaction : Interaction
Ph-ph interaction.
grid_points : array_like
Grid-point indices where imag-self-energeis are caclculated.
Grid-point indices where imag-self-energies are calculated.
dtype=int, shape=(grid_points,)
temperatures : array_like
Temperatures where imag-self-energies are calculated.
dtype=float, shape=(temperatures,)
sigmas : array_like, optional
sigmas : Sequence, optional
A set of sigmas. simgas=[None, ] means to use tetrahedron method,
otherwise smearing method with real positive value of sigma.
For example, sigmas=[None, 0.01, 0.03] is possible. Default is None,
@ -671,9 +674,7 @@ def get_imag_self_energy(
"""
if sigmas is None:
_sigmas = [
None,
]
_sigmas = [None]
else:
_sigmas = sigmas
@ -718,7 +719,7 @@ def get_imag_self_energy(
order="C",
)
detailed_gamma: List[Optional[np.ndarray]] = []
detailed_gamma = []
ise = ImagSelfEnergy(
interaction, with_detail=(write_gamma_detail or return_gamma_detail)
@ -771,10 +772,7 @@ def get_imag_self_energy(
log_level,
)
if return_gamma_detail:
return _frequency_points, gamma, detailed_gamma
else:
return _frequency_points, gamma
return _frequency_points, gamma, detailed_gamma
def _get_imag_self_energy_at_gp(

73
test/phonon3/test_imag_self_energy.py
View File

@ -25,7 +25,7 @@ def test_imag_self_energy_at_bands(si_pbesol: Phono3py):
[0.00382813, 0.0049497, 0.02727924, 0.01382784, 0.04133946, 0.02980282],
]
for i, grgp in enumerate((1, 103)):
_, gammas = si_pbesol.run_imag_self_energy(
ise_params = si_pbesol.run_imag_self_energy(
[
si_pbesol.grid.grg2bzg[grgp],
],
@ -35,7 +35,9 @@ def test_imag_self_energy_at_bands(si_pbesol: Phono3py):
frequency_points_at_bands=True,
)
# print(gammas.ravel())
np.testing.assert_allclose(gammas.ravel(), gammas_ref[i], rtol=0, atol=1e-2)
np.testing.assert_allclose(
ise_params.gammas.ravel(), gammas_ref[i], rtol=0, atol=1e-2
)
def test_imag_self_energy_at_bands_detailed(si_pbesol: Phono3py):
@ -47,7 +49,7 @@ def test_imag_self_energy_at_bands_detailed(si_pbesol: Phono3py):
"""
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
_, gammas, detailed_gammas = si_pbesol.run_imag_self_energy(
ise_params = si_pbesol.run_imag_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
@ -145,10 +147,14 @@ def test_imag_self_energy_at_bands_detailed(si_pbesol: Phono3py):
]
weights_103 = [2] * 364 + [1]
gammas_1_ref = gammas[:, :, 0].ravel()
gammas_103_ref = gammas[:, :, 1].ravel()
gammas_1 = np.dot(weights_1, detailed_gammas[0][0, 0].sum(axis=-1).sum(axis=-1))
gammas_103 = np.dot(weights_103, detailed_gammas[1][0, 0].sum(axis=-1).sum(axis=-1))
gammas_1_ref = ise_params.gammas[:, :, 0].ravel()
gammas_103_ref = ise_params.gammas[:, :, 1].ravel()
gammas_1 = np.dot(
weights_1, ise_params.detailed_gammas[0][0, 0].sum(axis=-1).sum(axis=-1)
)
gammas_103 = np.dot(
weights_103, ise_params.detailed_gammas[1][0, 0].sum(axis=-1).sum(axis=-1)
)
np.testing.assert_allclose(
gammas_1[:2].sum(), gammas_1_ref[:2].sum(), rtol=0, atol=1e-2
)
@ -477,7 +483,7 @@ def test_imag_self_energy_npoints(si_pbesol: Phono3py, with_given_freq_points: b
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
if with_given_freq_points:
fpoints, gammas = si_pbesol.run_imag_self_energy(
ise_params = si_pbesol.run_imag_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
@ -485,7 +491,7 @@ def test_imag_self_energy_npoints(si_pbesol: Phono3py, with_given_freq_points: b
frequency_points=ref_freq_points,
)
else:
fpoints, gammas = si_pbesol.run_imag_self_energy(
ise_params = si_pbesol.run_imag_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
@ -497,8 +503,12 @@ def test_imag_self_energy_npoints(si_pbesol: Phono3py, with_given_freq_points: b
# for line in gammas.reshape(-1, 10):
# print("[", ",".join([f"{val:.8f}" for val in line]), "],")
np.testing.assert_allclose(ref_gammas, gammas.reshape(-1, 10), rtol=0, atol=1e-2)
np.testing.assert_allclose(ref_freq_points, fpoints.ravel(), rtol=0, atol=1e-5)
np.testing.assert_allclose(
ref_gammas, ise_params.gammas.reshape(-1, 10), rtol=0, atol=1e-2
)
np.testing.assert_allclose(
ref_freq_points, ise_params.frequency_points.ravel(), rtol=0, atol=1e-5
)
def test_imag_self_energy_npoints_with_sigma(si_pbesol: Phono3py):
@ -820,7 +830,7 @@ def test_imag_self_energy_npoints_with_sigma(si_pbesol: Phono3py):
]
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
fpoints, gammas = si_pbesol.run_imag_self_energy(
ise_params = si_pbesol.run_imag_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
@ -832,8 +842,12 @@ def test_imag_self_energy_npoints_with_sigma(si_pbesol: Phono3py):
# for line in gammas.reshape(-1, 10):
# print("[", ",".join([f"{val:.8f}" for val in line]), "],")
np.testing.assert_allclose(ref_gammas, gammas.reshape(-1, 10), rtol=0, atol=1e-2)
np.testing.assert_allclose(ref_freq_points, fpoints.ravel(), rtol=0, atol=1e-5)
np.testing.assert_allclose(
ref_gammas, ise_params.gammas.reshape(-1, 10), rtol=0, atol=1e-2
)
np.testing.assert_allclose(
ref_freq_points, ise_params.frequency_points.ravel(), rtol=0, atol=1e-5
)
si_pbesol.sigmas = None
@ -884,7 +898,7 @@ def test_imag_self_energy_detailed(si_pbesol: Phono3py):
]
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
_, _, detailed_gammas = si_pbesol.run_imag_self_energy(
ise_params = si_pbesol.run_imag_self_energy(
si_pbesol.grid.grg2bzg[
[
1,
@ -898,13 +912,16 @@ def test_imag_self_energy_detailed(si_pbesol: Phono3py):
)
print(
",".join(
[f"{val:.8f}" for val in detailed_gammas[0][0, 0].sum(axis=(1, 2, 3, 4))]
[
f"{val:.8f}"
for val in ise_params.detailed_gammas[0][0, 0].sum(axis=(1, 2, 3, 4))
]
)
)
np.testing.assert_allclose(
ref_detailed_gamma,
detailed_gammas[0][0, 0].sum(axis=(1, 2, 3, 4)),
ise_params.detailed_gammas[0][0, 0].sum(axis=(1, 2, 3, 4)),
rtol=0,
atol=1e-2,
)
@ -1428,7 +1445,7 @@ def test_imag_self_energy_scat_classes(si_pbesol: Phono3py, scattering_class: in
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
_, gammas = si_pbesol.run_imag_self_energy(
ise_params = si_pbesol.run_imag_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
@ -1441,7 +1458,7 @@ def test_imag_self_energy_scat_classes(si_pbesol: Phono3py, scattering_class: in
np.testing.assert_allclose(
gammas_classes[scattering_class - 1],
gammas.ravel(),
ise_params.gammas.ravel(),
rtol=0,
atol=1e-2,
)
@ -1714,7 +1731,7 @@ def test_imag_self_energy_nacl_npoints(nacl_pbe: Phono3py):
nacl_pbe.mesh_numbers = [9, 9, 9]
nacl_pbe.init_phph_interaction()
fpoints, gammas = nacl_pbe.run_imag_self_energy(
ise_params = nacl_pbe.run_imag_self_energy(
nacl_pbe.grid.grg2bzg[[1, 103]],
[
300,
@ -1724,8 +1741,12 @@ def test_imag_self_energy_nacl_npoints(nacl_pbe: Phono3py):
# print(",".join([f"{val:.8f}" for val in gammas.ravel()]))
np.testing.assert_allclose(ref_gammas_nacl, gammas.ravel(), rtol=0, atol=2e-2)
np.testing.assert_allclose(ref_freq_points_nacl, fpoints.ravel(), rtol=0, atol=1e-5)
np.testing.assert_allclose(
ref_gammas_nacl, ise_params.gammas.ravel(), rtol=0, atol=2e-2
)
np.testing.assert_allclose(
ref_freq_points_nacl, ise_params.frequency_points.ravel(), rtol=0, atol=1e-5
)
def test_imag_self_energy_nacl_nac_npoints(nacl_pbe: Phono3py):
@ -1876,13 +1897,15 @@ def test_imag_self_energy_nacl_nac_npoints(nacl_pbe: Phono3py):
nacl_pbe.mesh_numbers = [9, 9, 9]
nacl_pbe.init_phph_interaction(nac_q_direction=[1, 0, 0])
fpoints, gammas = nacl_pbe.run_imag_self_energy(
ise_params = nacl_pbe.run_imag_self_energy(
[nacl_pbe.grid.gp_Gamma], [300], num_frequency_points=10
)
# print(",".join([f"{val:.8f}" for val in gammas.ravel()]))
np.testing.assert_allclose(
ref_freq_points_nacl_nac, fpoints.ravel(), rtol=0, atol=1e-5
ref_freq_points_nacl_nac, ise_params.frequency_points.ravel(), rtol=0, atol=1e-5
)
np.testing.assert_allclose(
ref_gammas_nacl_nac, ise_params.gammas.ravel(), rtol=0, atol=2e-2
)
np.testing.assert_allclose(ref_gammas_nacl_nac, gammas.ravel(), rtol=0, atol=2e-2)