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

Merge branch 'develop'

This commit is contained in:
Atsushi Togo 2025-07-22 16:15:29 +09:00
parent 7e6e2ce628 a12bb6becd
commit 8e7bc01229
39 changed files with 1060 additions and 595 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.pre-commit-config.yaml
View File

@ -11,7 +11,7 @@ repos:
exclude: ^example/AlN-LDA/
- repo: https://github.com/astral-sh/ruff-pre-commit
rev: v0.12.1
rev: v0.12.4
hooks:
- id: ruff
args: [ "--fix", "--show-fixes" ]

6
doc/changelog.md
View File

@ -2,6 +2,12 @@
# Change Log
## Jul-22-2025: Version 3.18.0
- Changed `Phono3py.run_imag_self_energy()` to return `ImagSelfEnergyValues`.
- Traditional force constants symmetrizer now applies translational and
permutation symmetries alternately 3 times in succession (previously once).
## 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)

4
doc/conf.py
View File

@ -58,9 +58,9 @@ copyright = "2015, Atsushi Togo"
# built documents.
#
# The short X.Y version.
version = "3.17"
version = "3.18"
# The full version, including alpha/beta/rc tags.
release = "3.17.1"
release = "3.18.0"
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.

63
example/AlN-rd/README.md Normal file
View File

@ -0,0 +1,63 @@
# AlN lattice thermal conductivity calculation from dataset for pypolymlp
## Computational setting of VASP calculations
For supercell forces and energies
- Supercell 4x4x2 of wurtzite unit cell
- Random directional displacements of 0.03 Angstrom
- PBE-sol
- 520 eV cutoff energy
- Gamma centered 2x2x2 kpoint mesh
- LREAL = .FALSE.
- ADDGRID = .TRUE.
For parameters of non-analytical term correction,
- PBE-sol
- 520 eV cutoff energy
- Gamma centered 7x7x4 kpoint mesh
- LEPSION = .TRUE.
- LREAL = .FALSE.
These data are stored in `phonopy_params_mp-661.yaml.xz`.
## Example of lattice thermal conductivity calculation
MLPs by pypolymlp are developed by
```bash
% phono3py-load phonopy_params_mp-661.yaml.xz --pypolymlp -v
```
Dataset with 180 supercells is used for training and 20 for the test. This
calculation will take 5-10 minutes depending on computer resource.
`pypolymlp.yaml` is made by this command.
Force constants are calculated by
```bash
% phono3py-load phonopy_params_mp-661.yaml.xz --pypolymlp --relax-atomic-positions -d
```
With the `--relax-atomic-positions` option, internal atomic positions in unit
cell are optimized by pypolymlp. The displacement-force dataset is stored in
`phono3py_mlp_eval_dataset.yaml`. Force constants are symmetried using symfc,
but the phono3py's traditional symmetrizer can be used with the option
`--fc-calculator traditional`. The symmetry constraints applied by this
traditional symmetrizer is weaker, but the calculation demands less memory
space.
Lattice thermal conductivity is calculated by
```bash
% phono3py-load phonopy_params_mp-661.yaml.xz --mesh 40 --br
```
Steps written above are performed in one-shot by
```bash
% phono3py-load phonopy_params_mp-661.yaml.xz --pypolymlp --relax-atomic-positions -d --mesh 40 --br
```
The lattice thermal conductivity calculated at 300 K will be around k_xx=252 and k_zz=232.

BIN
example/AlN-rd/phonopy_params_mp-661.yaml.xz Normal file
View File

Binary file not shown.

18
phono3py/__init__.py
View File

@ -34,8 +34,16 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from phono3py.api_isotope import Phono3pyIsotope # noqa F401
from phono3py.api_jointdos import Phono3pyJointDos # noqa F401
from phono3py.api_phono3py import Phono3py # noqa F401
from phono3py.cui.load import load # noqa F401
from phono3py.version import __version__ # noqa F401
from phono3py.api_isotope import Phono3pyIsotope
from phono3py.api_jointdos import Phono3pyJointDos
from phono3py.api_phono3py import Phono3py
from phono3py.cui.load import load
from phono3py.version import __version__
__all__ = [
"Phono3pyIsotope",
"Phono3pyJointDos",
"Phono3py",
"load",
"__version__",
]

3
phono3py/api_isotope.py
View File

@ -37,6 +37,7 @@
import numpy as np
from phono3py.other.isotope import Isotope
from phono3py.phonon.grid import BZGrid
class Phono3pyIsotope:
@ -81,7 +82,7 @@ class Phono3pyIsotope:
return self._iso.dynamical_matrix
@property
def grid(self):
def grid(self) -> BZGrid:
"""Return BZGrid class instance."""
return self._iso.bz_grid

7
phono3py/api_jointdos.py
View File

@ -34,7 +34,10 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import numpy as np
from phonopy.harmonic.dynamical_matrix import DynamicalMatrix
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import Primitive, Supercell
from phonopy.structure.symmetry import Symmetry
@ -116,7 +119,7 @@ class Phono3pyJointDos:
self.initialize(mesh)
@property
def grid(self):
def grid(self) -> BZGrid | None:
"""Return BZGrid class instance."""
return self._bz_grid
@ -293,7 +296,7 @@ class Phono3pyJointDos:
print('JDOS is written into "%s".' % filename)
@property
def dynamical_matrix(self):
def dynamical_matrix(self) -> DynamicalMatrix:
"""Return DynamicalMatrix class instance."""
return self._jdos.dynamical_matrix

356
phono3py/api_phono3py.py
View File

@ -36,12 +36,14 @@
from __future__ import annotations
import copy
import dataclasses
import os
import warnings
from collections.abc import Sequence
from typing import Literal, Optional, Union, cast
from typing import Literal, cast
import numpy as np
from numpy.typing import NDArray
from numpy.typing import ArrayLike, NDArray
from phonopy.harmonic.displacement import (
directions_to_displacement_dataset,
get_least_displacements,
@ -59,7 +61,11 @@ from phonopy.interface.mlp import PhonopyMLP
from phonopy.interface.pypolymlp import (
PypolymlpParams,
)
from phonopy.interface.symfc import SymfcFCSolver, symmetrize_by_projector
from phonopy.interface.symfc import (
SymfcFCSolver,
parse_symfc_options,
symmetrize_by_projector,
)
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import (
@ -75,7 +81,10 @@ from phonopy.structure.symmetry import Symmetry
from phono3py.conductivity.init_direct_solution import get_thermal_conductivity_LBTE
from phono3py.conductivity.init_rta import get_thermal_conductivity_RTA
from phono3py.interface.fc_calculator import FC3Solver
from phono3py.interface.fc_calculator import (
FC3Solver,
extract_fc2_fc3_calculators_options,
)
from phono3py.interface.phono3py_yaml import Phono3pyYaml
from phono3py.phonon.grid import BZGrid
from phono3py.phonon3.dataset import forces_in_dataset
@ -103,6 +112,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.
@ -147,19 +166,19 @@ class Phono3py:
def __init__(
self,
unitcell: PhonopyAtoms,
supercell_matrix=None,
primitive_matrix=None,
phonon_supercell_matrix=None,
cutoff_frequency=1e-4,
frequency_factor_to_THz=None,
is_symmetry=True,
is_mesh_symmetry=True,
use_grg=False,
SNF_coordinates="reciprocal",
supercell_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,
is_symmetry: bool = True,
is_mesh_symmetry: bool = True,
use_grg: bool = False,
SNF_coordinates: str = "reciprocal",
make_r0_average: bool = True,
symprec=1e-5,
calculator: Optional[str] = None,
log_level=0,
symprec: float = 1e-5,
calculator: str | None = None,
log_level: int = 0,
):
"""Init method.
@ -232,7 +251,7 @@ class Phono3py:
self._make_r0_average = make_r0_average
self._cutoff_frequency = cutoff_frequency
self._calculator: Optional[str] = calculator
self._calculator = calculator
self._log_level = log_level
# Create supercell and primitive cell
@ -240,8 +259,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(
@ -260,9 +278,7 @@ class Phono3py:
self._build_phonon_supercell()
self._build_phonon_primitive_cell()
self._sigmas = [
None,
]
self._sigmas = [None]
self._sigma_cutoff = None
# Grid
@ -277,18 +293,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
@ -305,9 +320,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
@ -326,7 +341,7 @@ class Phono3py:
return __version__
@property
def calculator(self) -> Optional[str]:
def calculator(self) -> str | None:
"""Return calculator interface name.
str
@ -424,7 +439,7 @@ class Phono3py:
self._sigmas.append(None)
@property
def sigma_cutoff(self) -> Optional[float]:
def sigma_cutoff(self) -> float | None:
"""Setter and getter of Smearing cutoff width.
This is given as a multiple of the standard deviation.
@ -441,7 +456,7 @@ class Phono3py:
self._sigma_cutoff = sigma_cutoff
@property
def nac_params(self) -> Optional[dict]:
def nac_params(self) -> dict | None:
"""Setter and getter of parameters for non-analytical term correction.
dict
@ -465,7 +480,7 @@ class Phono3py:
self._init_dynamical_matrix()
@property
def dynamical_matrix(self) -> Optional[DynamicalMatrix]:
def dynamical_matrix(self) -> DynamicalMatrix | None:
"""Return DynamicalMatrix instance.
This is not dynamical matrices but the instance of DynamicalMatrix
@ -582,10 +597,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'
@ -605,7 +620,7 @@ class Phono3py:
return self._frequency_factor_to_THz
@property
def dataset(self) -> Optional[dict]:
def dataset(self) -> dict | None:
"""Setter and getter of displacement-force dataset.
dict
@ -667,7 +682,7 @@ class Phono3py:
self._phonon_supercells_with_displacements = None
@property
def phonon_dataset(self) -> Optional[dict]:
def phonon_dataset(self) -> dict | None:
"""Setter and getter of displacement-force dataset for fc2.
dict
@ -712,7 +727,7 @@ class Phono3py:
self._phonon_supercells_with_displacements = None
@property
def mlp_dataset(self) -> Optional[dict]:
def mlp_dataset(self) -> dict | None:
"""Return displacement-force dataset.
The supercell matrix is equal to that of usual displacement-force
@ -728,7 +743,7 @@ class Phono3py:
self._mlp_dataset = mlp_dataset
@property
def phonon_mlp_dataset(self) -> Optional[dict]:
def phonon_mlp_dataset(self) -> dict | None:
"""Return phonon displacement-force dataset.
The phonon supercell matrix is equal to that of usual displacement-force
@ -803,7 +818,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
@ -813,6 +828,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
@ -834,7 +850,7 @@ class Phono3py:
return self._phonon_supercells_with_displacements
@property
def mesh_numbers(self):
def mesh_numbers(self) -> NDArray | None:
"""Setter and getter of sampling mesh numbers in reciprocal space."""
if self._bz_grid is None:
return None
@ -842,7 +858,7 @@ class Phono3py:
return self._bz_grid.D_diag
@mesh_numbers.setter
def mesh_numbers(self, mesh_numbers: Union[int, float, Sequence, NDArray]):
def mesh_numbers(self, mesh_numbers: float | ArrayLike):
self._set_mesh_numbers(mesh_numbers)
@property
@ -882,7 +898,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:
@ -1077,7 +1093,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):
@ -1091,12 +1109,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.
@ -1250,7 +1268,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,
@ -1395,12 +1413,12 @@ class Phono3py:
def generate_fc2_displacements(
self,
distance: Optional[float] = None,
is_plusminus: str = "auto",
distance: float | None = None,
is_plusminus: bool | Literal["auto"] = "auto",
is_diagonal: bool = False,
number_of_snapshots: Optional[Union[int, Literal["auto"]]] = None,
random_seed: Optional[int] = None,
max_distance: Optional[float] = None,
number_of_snapshots: int | Literal["auto"] | None = None,
random_seed: int | None = None,
max_distance: float | None = None,
):
"""Generate displacement dataset in phonon supercell for fc2.
@ -1541,9 +1559,13 @@ class Phono3py:
if fc_calculator == "traditional" or fc_calculator is None:
if symmetrize_fc3r:
self.symmetrize_fc3(
use_symfc_projector=use_symfc_projector and fc_calculator is None
fc3_calc_opts = extract_fc2_fc3_calculators_options(
fc_calculator_options, 3
)
if use_symfc_projector and fc_calculator is None:
self.symmetrize_fc3(use_symfc_projector=True, options=fc3_calc_opts)
else:
self.symmetrize_fc3(options=fc3_calc_opts)
elif fc_calculator == "symfc":
symfc_solver = cast(SymfcFCSolver, fc_solver.fc_solver)
fc3_nonzero_elems = symfc_solver.get_nonzero_atomic_indices_fc3()
@ -1566,36 +1588,76 @@ class Phono3py:
self._fc2 = fc2
if fc_calculator == "traditional" or fc_calculator is None:
if symmetrize_fc3r:
self.symmetrize_fc2(
use_symfc_projector=use_symfc_projector
and fc_calculator is None
fc2_calc_opts = extract_fc2_fc3_calculators_options(
fc_calculator_options, 2
)
if use_symfc_projector and fc_calculator is None:
self.symmetrize_fc2(
use_symfc_projector=True, options=fc2_calc_opts
)
else:
self.symmetrize_fc2(options=fc2_calc_opts)
def symmetrize_fc3(self, use_symfc_projector: bool = False):
"""Symmetrize fc3 by symfc projector or traditional approach."""
def symmetrize_fc3(
self,
use_symfc_projector: bool = False,
options: str | None = None,
):
"""Symmetrize fc3 by symfc projector or traditional approach.
Parameters
----------
use_symfc_projector : bool, optional
If True, the force constants are symmetrized by symfc projector
instead of traditional approach.
options : str or None, optional
For symfc projector:
"use_mkl=true" calls sparse_dot_mkl (required to install it).
For traditional symmetrization:
"level=N" applies translational and permutation symmetries
alternately N times in succession. Default is 3.
"""
if self._fc3 is None:
raise RuntimeError("fc3 is not set. Call produce_fc3 first.")
if use_symfc_projector:
if self._log_level:
print("Symmetrizing fc3 by symfc projector.", flush=True)
if options is None:
_options = None
else:
_options = parse_symfc_options(options, 3)
self._fc3 = symmetrize_by_projector(
self._supercell,
self._fc3,
3,
primitive=self._primitive,
options=_options,
log_level=self._log_level,
show_credit=True,
)
else:
level = 3
if options is not None:
for option in options.split(","):
if "level" in option.lower():
try:
level = int(option.split("=")[1].split()[0])
except ValueError:
pass
break
if self._log_level:
print("Symmetrizing fc3 by traditional approach.", flush=True)
if self._fc3.shape[0] == self._fc3.shape[1]:
set_translational_invariance_fc3(self._fc3)
set_permutation_symmetry_fc3(self._fc3)
else:
set_translational_invariance_compact_fc3(self._fc3, self._primitive)
set_permutation_symmetry_compact_fc3(self._fc3, self._primitive)
print(
f"Symmetrizing fc3 by traditional approach (N={level}).", flush=True
)
for _ in range(level):
if self._fc3.shape[0] == self._fc3.shape[1]:
set_translational_invariance_fc3(self._fc3)
set_permutation_symmetry_fc3(self._fc3)
else:
set_translational_invariance_compact_fc3(self._fc3, self._primitive)
set_permutation_symmetry_compact_fc3(self._fc3, self._primitive)
def produce_fc2(
self,
@ -1657,11 +1719,30 @@ class Phono3py:
if symmetrize_fc2 and (fc_calculator is None or fc_calculator == "traditional"):
self.symmetrize_fc2(
use_symfc_projector=use_symfc_projector and fc_calculator is None
use_symfc_projector=use_symfc_projector and fc_calculator is None,
options=fc_calculator_options,
)
def symmetrize_fc2(self, use_symfc_projector: bool = False):
"""Symmetrize fc2 by symfc projector or traditional approach."""
def symmetrize_fc2(
self,
use_symfc_projector: bool = False,
options: str | None = None,
):
"""Symmetrize fc2 by symfc projector or traditional approach.
Parameters
----------
use_symfc_projector : bool, optional
If True, the force constants are symmetrized by symfc projector
instead of traditional approach.
options : str or None, optional
For symfc projector:
"use_mkl=true" calls sparse_dot_mkl (required to install it).
For traditional symmetrization:
"level=N" applies translational and permutation symmetries
alternately N times in succession. Default is 3.
"""
if self._fc2 is None:
raise RuntimeError(
"fc2 is not set. Call produce_fc3 or produce_fc2 (phonon_fc2) first."
@ -1678,21 +1759,38 @@ class Phono3py:
if use_symfc_projector:
if self._log_level:
print("Symmetrizing fc2 by symfc projector.", flush=True)
if options is None:
_options = None
else:
_options = parse_symfc_options(options, 2)
self._fc2 = symmetrize_by_projector(
supercell,
self._fc2,
2,
primitive=primitive,
options=_options,
log_level=self._log_level,
show_credit=True,
)
else:
level = 3
if options is not None:
for option in options.split(","):
if "level" in option.lower():
try:
level = int(option.split("=")[1].split()[0])
except ValueError:
pass
break
if self._log_level:
print("Symmetrizing fc2 by traditional approach.", flush=True)
if self._fc2.shape[0] == self._fc2.shape[1]:
symmetrize_force_constants(self._fc2)
else:
symmetrize_compact_force_constants(self._fc2, primitive)
print(
f"Symmetrizing fc2 by traditional approach (N={level}).", flush=True
)
for _ in range(level):
if self._fc2.shape[0] == self._fc2.shape[1]:
symmetrize_force_constants(self._fc2)
else:
symmetrize_compact_force_constants(self._fc2, primitive)
def cutoff_fc3_by_zero(self, cutoff_distance, fc3=None):
"""Set zero to fc3 elements out of cutoff distance.
@ -1753,7 +1851,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.
@ -1818,7 +1916,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,
@ -1836,25 +1933,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,
@ -2039,16 +2145,16 @@ class Phono3py:
def run_thermal_conductivity(
self,
is_LBTE: bool = False,
temperatures: Optional[Sequence] = None,
temperatures: Sequence | None = None,
is_isotope: bool = False,
mass_variances: Optional[Sequence] = None,
grid_points: Optional[Sequence[int]] = None,
boundary_mfp: Optional[float] = None, # in micrometer
mass_variances: Sequence | None = None,
grid_points: ArrayLike | None = None,
boundary_mfp: float | None = None, # in micrometer
solve_collective_phonon: bool = False,
use_ave_pp: bool = False,
is_reducible_collision_matrix: bool = False,
is_kappa_star: bool = True,
gv_delta_q: Optional[float] = None, # for group velocity
gv_delta_q: float | None = None, # for group velocity
is_full_pp: bool = False,
pinv_cutoff: float = 1.0e-8, # for pseudo-inversion of collision matrix
pinv_method: int = 0, # for pseudo-inversion of collision matrix
@ -2056,18 +2162,18 @@ class Phono3py:
write_gamma: bool = False,
read_gamma: bool = False,
is_N_U: bool = False,
conductivity_type: Optional[str] = None,
conductivity_type: str | None = None,
write_kappa: bool = False,
write_gamma_detail: bool = False,
write_collision: bool = False,
read_collision: bool = False,
read_collision: str | Sequence | None = None,
write_pp: bool = False,
read_pp: bool = False,
write_LBTE_solution: bool = False,
compression: str = "gzip",
input_filename: Optional[str] = None,
output_filename: Optional[str] = None,
log_level: Optional[int] = None,
input_filename: str | None = None,
output_filename: str | None = None,
log_level: int | None = None,
):
"""Run thermal conductivity calculation.
@ -2174,9 +2280,9 @@ class Phono3py:
is written into a file. With multiple `sigmas` specified,
respective files are created. Be careful that this file can be
huge.
read_collision : bool, optional, default is False
Direct solution only (`is_LBTE=True`). With True, collision matrix
is read from a file.
read_collision : str | Sequence, optional, default is None.
Direct solution only (`is_LBTE=True`). With specified, collision
matrix is read from a file.
write_pp : bool, optional, default is False
With True, phonon-phonon interaction strength is written into
files at each grid point. This option assumes single value is in
@ -2294,7 +2400,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.
@ -2348,14 +2456,14 @@ class Phono3py:
test_size=test_size,
)
def save_mlp(self, filename: Optional[str] = None):
def save_mlp(self, filename: str | None = None):
"""Save machine learning potential."""
if self._mlp is None:
raise RuntimeError("MLP is not developed yet.")
self._mlp.save(filename=filename)
def load_mlp(self, filename: Optional[str] = None):
def load_mlp(self, filename: str | os.PathLike | None = None):
"""Load machine learning potential."""
self._mlp = PhonopyMLP(log_level=self._log_level)
self._mlp.load(filename=filename)
@ -2385,7 +2493,7 @@ class Phono3py:
def develop_phonon_mlp(
self,
params: Optional[Union[PypolymlpParams, dict, str]] = None,
params: PypolymlpParams | dict | str | None = None,
test_size: float = 0.1,
):
"""Develop MLP for fc2.
@ -2412,14 +2520,14 @@ class Phono3py:
test_size=test_size,
)
def save_phonon_mlp(self, filename: Optional[str] = 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)
def load_phonon_mlp(self, filename: Optional[str] = None):
def load_phonon_mlp(self, filename: str | None = None):
"""Load machine learning potential."""
self._phonon_mlp = PhonopyMLP(log_level=self._log_level)
self._phonon_mlp.load(filename=filename)
@ -2437,8 +2545,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(
@ -2575,15 +2685,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"]
@ -2621,7 +2738,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)
@ -2630,7 +2749,7 @@ class Phono3py:
def _set_mesh_numbers(
self,
mesh: Union[int, float, Sequence, NDArray],
mesh: float | ArrayLike,
):
# initialization related to mesh
self._interaction = None
@ -2661,7 +2780,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):
@ -2717,6 +2837,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"]:
@ -2800,13 +2923,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,
@ -2815,7 +2936,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):

10
phono3py/conductivity/base.py
View File

@ -41,7 +41,7 @@ from abc import ABC, abstractmethod
from typing import Optional
import numpy as np
from numpy.typing import NDArray
from numpy.typing import ArrayLike, NDArray
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.phonon.thermal_properties import mode_cv
from phonopy.physical_units import get_physical_units
@ -419,12 +419,12 @@ class ConductivityBase(ABC):
def __init__(
self,
interaction: Interaction,
grid_points: np.ndarray | None = None,
temperatures: list | np.ndarray | None = None,
sigmas: list | np.ndarray | None = None,
grid_points: ArrayLike | None = None,
temperatures: ArrayLike | None = None,
sigmas: ArrayLike | None = None,
sigma_cutoff: float | None = None,
is_isotope=False,
mass_variances: list | np.ndarray | None = None,
mass_variances: ArrayLike | None = None,
boundary_mfp: float | None = None,
is_kappa_star: bool = True,
is_full_pp: bool = False,

19
phono3py/conductivity/direct_solution.py
View File

@ -36,9 +36,10 @@
from __future__ import annotations
from typing import Optional, Union
import os
import numpy as np
from numpy.typing import ArrayLike
from phono3py.conductivity.base import ConductivityComponents
from phono3py.conductivity.direct_solution_base import (
@ -54,20 +55,20 @@ class ConductivityLBTE(ConductivityLBTEBase):
def __init__(
self,
interaction: Interaction,
grid_points: Optional[np.ndarray] = None,
temperatures: Optional[Union[list, np.ndarray]] = None,
sigmas: Optional[Union[list, np.ndarray]] = None,
sigma_cutoff: Optional[float] = None,
grid_points: ArrayLike | None = None,
temperatures: ArrayLike | None = None,
sigmas: ArrayLike | None = None,
sigma_cutoff: float | None = None,
is_isotope: bool = False,
mass_variances: Optional[Union[list, np.ndarray]] = None,
boundary_mfp: Optional[float] = None, # in micrometer
mass_variances: ArrayLike | None = None,
boundary_mfp: float | None = None, # in micrometer
solve_collective_phonon: bool = False,
is_reducible_collision_matrix: bool = False,
is_kappa_star: bool = True,
gv_delta_q: Optional[float] = None,
gv_delta_q: float | None = None,
is_full_pp: bool = False,
read_pp: bool = False,
pp_filename: Optional[float] = None,
pp_filename: str | os.PathLike | None = None,
pinv_cutoff: float = 1.0e-8,
pinv_solver: int = 0,
pinv_method: int = 0,

18
phono3py/conductivity/direct_solution_base.py
View File

@ -36,12 +36,14 @@
from __future__ import annotations
import os
import sys
import time
from abc import abstractmethod
from typing import Optional, Union
from typing import Optional
import numpy as np
from numpy.typing import ArrayLike
from phonopy.phonon.degeneracy import degenerate_sets
from phonopy.physical_units import get_physical_units
@ -63,19 +65,19 @@ class ConductivityLBTEBase(ConductivityBase):
def __init__(
self,
interaction: Interaction,
grid_points: Optional[np.ndarray] = None,
temperatures: Optional[Union[list, np.ndarray]] = None,
sigmas: Optional[Union[list, np.ndarray]] = None,
sigma_cutoff: Optional[float] = None,
grid_points: ArrayLike | None = None,
temperatures: ArrayLike | None = None,
sigmas: ArrayLike | None = None,
sigma_cutoff: float | None = None,
is_isotope: bool = False,
mass_variances: Optional[Union[list, np.ndarray]] = None,
boundary_mfp: Optional[float] = None, # in micrometer
mass_variances: ArrayLike | None = None,
boundary_mfp: float | None = None, # in micrometer
solve_collective_phonon: bool = False,
is_reducible_collision_matrix: bool = False,
is_kappa_star: bool = True,
is_full_pp: bool = False,
read_pp: bool = False,
pp_filename: Optional[float] = None,
pp_filename: str | os.PathLike | None = None,
pinv_cutoff: float = 1.0e-8,
pinv_solver: int = 0,
pinv_method: int = 0,

66
phono3py/conductivity/init_direct_solution.py
View File

@ -34,8 +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 os
import sys
from typing import Optional, Union
from collections.abc import Sequence
from typing import Union
from numpy.typing import ArrayLike
from phono3py.conductivity.base import get_unit_to_WmK
from phono3py.conductivity.direct_solution import ConductivityLBTE
@ -59,32 +65,32 @@ cond_LBTE_type = Union[ConductivityLBTE, ConductivityWignerLBTE]
def get_thermal_conductivity_LBTE(
interaction: Interaction,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
is_isotope=False,
mass_variances=None,
grid_points=None,
boundary_mfp=None, # in micrometer
solve_collective_phonon=False,
is_reducible_collision_matrix=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
conductivity_type=None,
pinv_cutoff=1.0e-8,
pinv_solver=0, # default: dsyev in lapacke
pinv_method=0, # default: abs(eig) < cutoff
write_collision=False,
read_collision=False,
write_kappa=False,
write_pp=False,
read_pp=False,
write_LBTE_solution=False,
compression="gzip",
input_filename=None,
output_filename=None,
log_level=0,
temperatures: Sequence | None = None,
sigmas: Sequence | None = None,
sigma_cutoff: float | None = None,
is_isotope: bool = False,
mass_variances: Sequence | None = None,
grid_points: ArrayLike | None = None,
boundary_mfp: float | None = None, # in micrometer
solve_collective_phonon: bool = False,
is_reducible_collision_matrix: bool = False,
is_kappa_star: bool = True,
gv_delta_q: float | None = None,
is_full_pp: bool = False,
conductivity_type: str | None = None,
pinv_cutoff: float = 1.0e-8,
pinv_solver: int = 0, # default: dsyev in lapacke
pinv_method: int = 0, # default: abs(eig) < cutoff
write_collision: bool = False,
read_collision: str | Sequence | None = None,
write_kappa: bool = False,
write_pp: bool = False,
read_pp: bool = False,
write_LBTE_solution: bool = False,
compression: str = "gzip",
input_filename: str | os.PathLike | None = None,
output_filename: str | os.PathLike | None = None,
log_level: int = 0,
):
"""Calculate lattice thermal conductivity by direct solution."""
if temperatures is None:
@ -292,9 +298,9 @@ class ConductivityLBTEWriter:
volume: float,
is_reducible_collision_matrix: bool = False,
write_LBTE_solution: bool = False,
pinv_solver: Optional[int] = None,
pinv_solver: int | None = None,
compression: str = "gzip",
filename: Optional[str] = None,
filename: str | os.PathLike | None = None,
log_level: int = 0,
):
"""Write kappa related properties into a hdf5 file."""
@ -473,7 +479,7 @@ class ConductivityLBTEWriter:
def _set_collision_from_file(
lbte: ConductivityLBTEBase,
indices="all",
indices: str | Sequence | None = "all",
is_reducible_collision_matrix=False,
filename=None,
log_level=0,

2
phono3py/conductivity/wigner_direct_solution.py
View File

@ -34,6 +34,8 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import numpy as np
from phonopy.physical_units import get_physical_units

128
phono3py/cui/create_force_constants.py
View File

@ -77,7 +77,7 @@ def parse_forces(
without writing calculator name in it.
"""
filename_read_from: str | None = None
filename_read_from = None
dataset = None
if phono3py.phonon_supercell is None or fc_type == "fc3":
@ -170,7 +170,9 @@ def _read_FORCES_FC3_or_FC2(
print(f'{n_disp} snapshots were found in "{filename}".')
return _dataset
# Type-1
# Try reading type-1 dataset
if dataset is None:
raise RuntimeError("Type-1 displacement dataset is not given.")
if fc_type == "fc3":
parse_FORCES_FC3(dataset, filename)
else:
@ -183,18 +185,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,
):
@ -219,12 +212,12 @@ def run_pypolymlp_to_compute_forces(
".xz",
".gz",
".bz2",
"lzma",
".lzma",
]:
continue
if log_level:
print(f'Load MLPs from "{mlp_filename}".')
ph3py.load_mlp(mlp_filename)
print(f'Load MLPs from "{_mlp_filename}".')
ph3py.load_mlp(_mlp_filename)
mlp_loaded = True
if log_level and mlp_filename == "phono3py.pmlp":
print(f'Loading MLPs from "{_mlp_filename}" is obsolete.')
@ -253,58 +246,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:

1
phono3py/cui/kaccum_script.py
View File

@ -229,6 +229,7 @@ def _read_files(args: argparse.Namespace) -> tuple[h5py.File, PhonopyAtoms | Non
cell_info = collect_cell_info(
supercell_matrix=np.eye(3, dtype=int),
phonopy_yaml_cls=Phono3pyYaml,
load_phonopy_yaml=True,
)
cell_filename = cell_info.optional_structure_info[0]
print(f'# Crystal structure was read from "{cell_filename}".')

16
phono3py/cui/load.py
View File

@ -51,12 +51,13 @@ 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 (
extract_fc2_fc3_calculators,
extract_fc2_fc3_calculators_options,
update_cutoff_fc_calculator_options,
)
from phono3py.interface.phono3py_yaml import Phono3pyYaml
@ -292,6 +293,7 @@ def load(
elif phono3py_yaml is not None:
ph3py_yaml = Phono3pyYaml()
ph3py_yaml.read(phono3py_yaml)
assert ph3py_yaml.unitcell is not None
cell = ph3py_yaml.unitcell.copy()
_calculator = ph3py_yaml.calculator
smat = ph3py_yaml.supercell_matrix
@ -366,14 +368,13 @@ def load(
)
if use_pypolymlp and ph3py.fc3 is None and forces_in_dataset(ph3py.dataset):
assert ph3py.dataset is not None
ph3py.mlp_dataset = ph3py.dataset
ph3py.dataset = None
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,
@ -432,12 +433,11 @@ def compute_force_constants_from_datasets(
"""
fc3_calculator = extract_fc2_fc3_calculators(fc_calculator, 3)
fc2_calculator = extract_fc2_fc3_calculators(fc_calculator, 2)
fc3_calc_opts = extract_fc2_fc3_calculators(fc_calculator_options, 3)
fc3_calc_opts = extract_fc2_fc3_calculators_options(fc_calculator_options, 3)
fc3_calc_opts = update_cutoff_fc_calculator_options(
fc3_calc_opts, cutoff_pair_distance
)
fc2_calc_opts = extract_fc2_fc3_calculators(fc_calculator_options, 2)
exist_fc2 = ph3py.fc2 is not None
fc2_calc_opts = extract_fc2_fc3_calculators_options(fc_calculator_options, 2)
if ph3py.fc3 is None and forces_in_dataset(ph3py.dataset):
ph3py.produce_fc3(
symmetrize_fc3r=symmetrize_fc,
@ -447,7 +447,7 @@ def compute_force_constants_from_datasets(
use_symfc_projector=load_phono3py_yaml,
)
if not exist_fc2:
if ph3py.fc2 is None or fc3_calculator != fc2_calculator:
if (
ph3py.phonon_supercell_matrix is None and forces_in_dataset(ph3py.dataset)
) or (

55
phono3py/cui/phono3py_argparse.py
View File

@ -34,6 +34,8 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import sys
from phonopy.cui.phonopy_argparse import fix_deprecated_option_names
@ -56,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(
@ -178,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(
@ -186,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:
@ -278,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(
@ -336,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",
)
@ -344,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 "
@ -381,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(
@ -404,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"
),
@ -427,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(
@ -544,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:
@ -593,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",
@ -739,6 +727,13 @@ def get_parser(load_phono3py_yaml: bool = False):
default=None,
help="Solve reducible collision matrix",
)
parser.add_argument(
"--relax-atomic-positions",
dest="relax_atomic_positions",
action="store_true",
default=None,
help="Relax atomic positions using polynomial MLPs",
)
parser.add_argument(
"--rse",
dest="is_real_self_energy",

379
phono3py/cui/phono3py_script.py
View File

@ -40,12 +40,11 @@ import argparse
import os
import pathlib
import sys
import warnings
from collections.abc import Sequence
from typing import Optional, cast
from typing import cast
import numpy as np
from numpy.typing import NDArray
from numpy.typing import ArrayLike, NDArray
from phonopy.api_phonopy import Phonopy
from phonopy.cui.phonopy_argparse import show_deprecated_option_warnings
from phonopy.cui.phonopy_script import (
@ -58,17 +57,27 @@ from phonopy.cui.phonopy_script import (
store_nac_params,
)
from phonopy.cui.settings import PhonopySettings
from phonopy.exception import CellNotFoundError, ForceCalculatorRequiredError
from phonopy.exception import (
CellNotFoundError,
ForceCalculatorRequiredError,
PypolymlpRelaxationError,
)
from phonopy.file_IO import is_file_phonopy_yaml
from phonopy.harmonic.dynamical_matrix import DynamicalMatrixGL
from phonopy.harmonic.force_constants import show_drift_force_constants
from phonopy.interface.calculator import get_calculator_physical_units
from phonopy.interface.pypolymlp import get_change_in_positions, relax_atomic_positions
from phonopy.interface.symfc import estimate_symfc_cutoff_from_memsize
from phonopy.phonon.band_structure import get_band_qpoints
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
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,
@ -192,7 +201,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:
@ -215,6 +223,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
@ -235,7 +245,7 @@ def _start_phono3py(**argparse_control) -> tuple[argparse.Namespace, int]:
# Title
if log_level:
print_phono3py()
import phono3py._phono3py as phono3c
import phono3py._phono3py as phono3c # type: ignore[import]
max_threads = phono3c.omp_max_threads()
if max_threads > 0:
@ -313,33 +323,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
@ -454,7 +437,7 @@ def _check_supercell_in_yaml(
sys.exit(1)
def _init_phono3py(
def _init_phono3py_with_cell_info(
settings: Phono3pySettings,
cell_info: Phono3pyCellInfoResult,
interface_mode: str | None,
@ -462,28 +445,54 @@ def _init_phono3py(
log_level: int,
) -> tuple[Phono3py, dict]:
"""Initialize phono3py and update settings by default values."""
physical_units = get_calculator_physical_units(interface_mode)
distance_to_A = physical_units["distance_to_A"]
# Change unit of lattice parameters to angstrom
unitcell = cell_info.unitcell.copy()
if distance_to_A is not None:
lattice = unitcell.cell
lattice *= distance_to_A
unitcell.cell = lattice
# updated_settings keys
# ('sigmas', 'temperature_points', 'temperatures',
# 'frequency_factor_to_THz', 'num_frequency_points',
# 'frequency_step', 'frequency_scale_factor',
# 'cutoff_frequency')
phono3py, updated_settings, distance_to_A = _init_phono3py(
settings,
cell_info.unitcell.copy(),
supercell_matrix=cell_info.supercell_matrix,
primitive_matrix=cell_info.primitive_matrix,
phonon_supercell_matrix=cell_info.phonon_supercell_matrix,
interface_mode=interface_mode,
symprec=symprec,
log_level=log_level,
)
_check_supercell_in_yaml(cell_info, phono3py, distance_to_A, log_level)
return phono3py, updated_settings
def _init_phono3py(
settings: Phono3pySettings,
unitcell: PhonopyAtoms,
supercell_matrix: ArrayLike | None = None,
primitive_matrix: ArrayLike | None = None,
phonon_supercell_matrix: ArrayLike | None = None,
interface_mode: str | None = None,
symprec: float = 1e-5,
log_level: int = 0,
) -> tuple[Phono3py, dict, float | None]:
"""Initialize phono3py and update settings by default values."""
if interface_mode is not None:
physical_units = get_calculator_physical_units(interface_mode)
distance_to_A = physical_units["distance_to_A"]
assert distance_to_A is not None
lattice = unitcell.cell
lattice *= distance_to_A
unitcell.cell = lattice
else:
distance_to_A = None
updated_settings = _get_default_values(settings)
phono3py = Phono3py(
unitcell,
cell_info.supercell_matrix,
primitive_matrix=cell_info.primitive_matrix,
phonon_supercell_matrix=cell_info.phonon_supercell_matrix,
supercell_matrix,
primitive_matrix=primitive_matrix,
phonon_supercell_matrix=phonon_supercell_matrix,
cutoff_frequency=updated_settings["cutoff_frequency"],
frequency_factor_to_THz=updated_settings["frequency_factor_to_THz"],
is_symmetry=settings.is_symmetry,
@ -498,12 +507,12 @@ def _init_phono3py(
phono3py.sigmas = updated_settings["sigmas"]
phono3py.sigma_cutoff = settings.sigma_cutoff_width
_check_supercell_in_yaml(cell_info, phono3py, distance_to_A, log_level)
return phono3py, updated_settings
return phono3py, updated_settings, distance_to_A
def _settings_to_grid_points(settings: Phono3pySettings, bz_grid: BZGrid):
def _settings_to_grid_points(
settings: Phono3pySettings, bz_grid: BZGrid
) -> ArrayLike | None:
"""Read or set grid point indices."""
if settings.grid_addresses is not None:
grid_points = _grid_addresses_to_grid_points(settings.grid_addresses, bz_grid)
@ -514,7 +523,7 @@ def _settings_to_grid_points(settings: Phono3pySettings, bz_grid: BZGrid):
return grid_points
def _grid_addresses_to_grid_points(grid_addresses: NDArray, bz_grid: BZGrid):
def _grid_addresses_to_grid_points(grid_addresses: NDArray, bz_grid: BZGrid) -> NDArray:
"""Return grid point indices from grid addresses."""
grid_points = [
get_grid_point_from_address(ga, bz_grid.D_diag) for ga in grid_addresses
@ -527,9 +536,10 @@ def _create_supercells_with_displacements(
cell_info: Phono3pyCellInfoResult,
confs_dict: dict,
unitcell_filename: str,
interface_mode: Optional[str],
symprec: float,
log_level: int,
interface_mode: str | None,
symprec: float = 1e-5,
output_yaml_filename: str | None = None,
log_level: int = 0,
):
"""Create supercells and write displacements."""
if (
@ -537,7 +547,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,
@ -547,7 +557,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,
@ -556,23 +566,123 @@ 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,
) -> Phono3py:
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,
)
_ph3py = ph3py
if settings.relax_atomic_positions:
if log_level:
print("Relaxing atomic positions using polynomial MLPs...")
try:
assert ph3py.mlp is not None
relaxed_unitcell = relax_atomic_positions(
ph3py.unitcell,
ph3py.mlp.mlp,
verbose=log_level > 1,
)
except (PypolymlpRelaxationError, ValueError) as e:
# ValueError can come from pypolymlp directly.
print_error_message(str(e))
if log_level:
print_error()
sys.exit(1)
if log_level:
if relaxed_unitcell is None:
print("No relaxation was performed due to symmetry constraints.")
else:
get_change_in_positions(
relaxed_unitcell, ph3py.unitcell, verbose=log_level > 0
)
print("-" * 76)
if relaxed_unitcell is not None:
if log_level:
print(
"Relaxed crystal structures will be stored in "
f'"{output_yaml_filename}".'
)
_ph3py, _, _ = _init_phono3py(
settings,
relaxed_unitcell.copy(),
supercell_matrix=ph3py.supercell_matrix,
primitive_matrix=ph3py.primitive_matrix,
phonon_supercell_matrix=ph3py.phonon_supercell_matrix,
symprec=ph3py.symmetry.tolerance,
log_level=log_level,
)
if ph3py.mesh_numbers is not None:
assert settings.mesh_numbers is not None
_ph3py.mesh_numbers = settings.mesh_numbers
_ph3py.nac_params = ph3py.nac_params
_ph3py.mlp = ph3py.mlp
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
)
return _ph3py
def _produce_force_constants(
ph3py: Phono3py,
settings: Phono3pySettings,
@ -622,7 +732,8 @@ def _produce_force_constants(
except ForceCalculatorRequiredError as e:
if load_phono3py_yaml:
if log_level:
print("Symfc will be used to handle general (or random) displacements.")
print(str(e))
print("Try symfc to handle general (or random) displacements.")
else:
print_error_message(str(e))
if log_level:
@ -690,7 +801,9 @@ def _produce_force_constants(
def _run_gruneisen_then_exit(
phono3py: Phono3py, settings: Phono3pySettings, output_filename: str, log_level: int
phono3py: Phono3py,
settings: Phono3pySettings,
log_level: int,
):
"""Run mode Grueneisen parameter calculation from fc3."""
if (
@ -703,6 +816,8 @@ def _run_gruneisen_then_exit(
print_error()
sys.exit(1)
assert phono3py.fc2 is not None
assert phono3py.fc3 is not None
if len(phono3py.fc2) != len(phono3py.fc3):
print("Supercells used for fc2 and fc3 have to be same.")
if log_level:
@ -734,7 +849,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,
)
@ -747,7 +861,6 @@ def _run_jdos_then_exit(
phono3py: Phono3py,
settings: Phono3pySettings,
updated_settings: dict,
output_filename: str | None,
log_level: int,
):
"""Run joint-DOS calculation."""
@ -769,15 +882,15 @@ 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,
)
if log_level > 0:
dm = joint_dos.dynamical_matrix
if dm.is_nac() and dm.nac_method == "gonze":
dm.show_Gonze_nac_message()
if isinstance(dm, DynamicalMatrixGL):
dm.show_nac_message()
assert joint_dos.grid is not None
grid_points = _settings_to_grid_points(settings, joint_dos.grid)
joint_dos.run(grid_points, write_jdos=True)
@ -819,8 +932,8 @@ def _run_isotope_then_exit(
)
if log_level > 0:
dm = iso.dynamical_matrix
if dm.is_nac() and dm.nac_method == "gonze":
dm.show_Gonze_nac_message()
if isinstance(dm, DynamicalMatrixGL):
dm.show_nac_message()
grid_points = _settings_to_grid_points(settings, iso.grid)
iso.run(grid_points)
@ -834,14 +947,13 @@ 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."""
if log_level:
print("Generating grid system ... ", end="", flush=True)
phono3py.mesh_numbers = settings.mesh_numbers
assert phono3py.grid is not None
assert phono3py.mesh_numbers is not None
bz_grid = phono3py.grid
if log_level:
if bz_grid.grid_matrix is None:
@ -872,7 +984,7 @@ def _init_phph_interaction(
if log_level:
print("-" * 27 + " Phonon calculations " + "-" * 28)
dm = phono3py.dynamical_matrix
if dm.is_nac() and dm.nac_method == "gonze":
if isinstance(dm, DynamicalMatrixGL):
dm.show_nac_message()
print("Running harmonic phonon calculations...")
sys.stdout.flush()
@ -891,7 +1003,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:
@ -903,9 +1014,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:
@ -971,12 +1080,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(
@ -990,7 +1097,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()
@ -1045,6 +1152,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 #
@ -1056,8 +1164,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,
)
#######################
@ -1068,7 +1177,7 @@ def main(**argparse_control):
# 'frequency_factor_to_THz', 'num_frequency_points',
# 'frequency_step', 'frequency_scale_factor',
# 'cutoff_frequency')
ph3py, updated_settings = _init_phono3py(
ph3py, updated_settings = _init_phono3py_with_cell_info(
settings, cell_info, interface_mode, symprec, log_level
)
@ -1081,8 +1190,6 @@ def main(**argparse_control):
run_mode,
ph3py,
unitcell_filename,
input_filename,
output_filename,
interface_mode,
)
@ -1131,6 +1238,7 @@ def main(**argparse_control):
"show_triplets_info",
)
run_modes_with_gp = ("imag_self_energy", "real_self_energy", "jdos", "isotope")
if settings.mesh_numbers is None and run_mode in run_modes_with_mesh:
print("")
print("Mesh numbers have to be specified.")
@ -1151,11 +1259,20 @@ def main(**argparse_control):
print_error()
sys.exit(1)
####################
# Set mesh numbers #
####################
if run_mode in run_modes_with_mesh:
assert settings.mesh_numbers is not None
# jdos and isotope modes need to set mesh numbers differently.
if run_mode not in ("jdos", "isotope"):
ph3py.mesh_numbers = settings.mesh_numbers
#########################################################
# Write ir-grid points and grid addresses and then exit #
#########################################################
if run_mode == "write_grid_info":
ph3py.mesh_numbers = settings.mesh_numbers
assert ph3py.grid is not None
write_grid_points(
ph3py.primitive,
ph3py.grid,
@ -1175,7 +1292,7 @@ def main(**argparse_control):
# Show reduced number of triplets at grid points and then exit #
################################################################
if run_mode == "show_triplets_info":
ph3py.mesh_numbers = settings.mesh_numbers
assert ph3py.grid is not None
grid_points = _settings_to_grid_points(settings, ph3py.grid)
show_num_triplets(
ph3py.primitive,
@ -1228,48 +1345,23 @@ 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(
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,
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
if ph3py.fc3 is None or (
ph3py.fc2 is None and ph3py.phonon_supercell_matrix is None
):
# Note that ph3py is replaced when relax_atomic_positions is True.
ph3py = _run_pypolymlp(
ph3py,
settings,
confs_dict,
output_yaml_filename=output_yaml_filename,
log_level=log_level,
)
else:
if log_level:
print(
"Pypolymlp was not developed or used because fc2 and fc3 "
"are available."
)
###########################
# Produce force constants #
@ -1280,21 +1372,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 #
@ -1324,13 +1414,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,
)
@ -1338,6 +1426,7 @@ def main(**argparse_control):
# Run imaginary part of self energy of bubble diagram #
#######################################################
if run_mode == "imag_self_energy":
assert ph3py.grid is not None
ph3py.run_imag_self_energy(
_settings_to_grid_points(settings, ph3py.grid),
updated_settings["temperature_points"],
@ -1347,13 +1436,13 @@ 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,
)
#####################################################
# Run frequency shift calculation of bubble diagram #
#####################################################
elif run_mode == "real_self_energy":
assert ph3py.grid is not None
ph3py.run_real_self_energy(
_settings_to_grid_points(settings, ph3py.grid),
updated_settings["temperature_points"],
@ -1361,13 +1450,13 @@ def main(**argparse_control):
num_frequency_points=updated_settings["num_frequency_points"],
write_txt=True,
write_hdf5=True,
output_filename=output_filename,
)
#######################################################
# Run spectral function calculation of bubble diagram #
#######################################################
elif run_mode == "spectral_function":
assert ph3py.grid is not None
ph3py.run_spectral_function(
_settings_to_grid_points(settings, ph3py.grid),
updated_settings["temperature_points"],
@ -1376,13 +1465,13 @@ def main(**argparse_control):
num_points_in_batch=updated_settings["num_points_in_batch"],
write_txt=True,
write_hdf5=True,
output_filename=output_filename,
)
####################################
# Run lattice thermal conductivity #
####################################
elif run_mode == "conductivity-RTA" or run_mode == "conductivity-LBTE":
assert ph3py.grid is not None
grid_points = _settings_to_grid_points(settings, ph3py.grid)
ph3py.run_thermal_conductivity(
is_LBTE=settings.is_lbte,
@ -1412,8 +1501,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:

3
phono3py/cui/settings.py
View File

@ -38,6 +38,7 @@ from __future__ import annotations
import argparse
import os
from typing import Literal
import numpy as np
from phonopy.cui.settings import ConfParser, Settings, fracval
@ -81,7 +82,7 @@ class Phono3pySettings(Settings):
self.is_symmetrize_fc3_q = False
self.is_symmetrize_fc3_r = False
self.is_tetrahedron_method = False
self.lapack_zheev_uplo = "L"
self.lapack_zheev_uplo: Literal["L", "U"] = "L"
self.mass_variances = None
self.max_freepath = None
self.num_points_in_batch = None

41
phono3py/cui/show_log.py
View File

@ -34,9 +34,10 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import sys
from collections.abc import Sequence
from typing import Optional, Union
import numpy as np
from phonopy.structure.cells import print_cell
@ -46,13 +47,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 = (
@ -64,12 +63,13 @@ def show_general_settings(
phonon_supercell_matrix = phono3py.phonon_supercell_matrix
print("-" * 29 + " General settings " + "-" * 29)
if run_mode:
if settings.use_pypolymlp:
if settings.create_displacements or settings.random_displacements is not None:
print(f"Run mode: pypolymlp + {run_mode}")
else:
print("Run mode: pypolymlp")
else:
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.hdf5_compression:
print(f"HDF5 data compression filter: {settings.hdf5_compression}")
if interface_mode:
@ -82,19 +82,18 @@ def show_general_settings(
print_supercell_matrix(supercell_matrix, phonon_supercell_matrix)
if is_primitive_axes_auto:
print("Primitive matrix (Auto):")
for v in primitive_matrix:
print(f" {v}")
elif primitive_matrix is not None:
print("Primitive matrix:")
if primitive_matrix is not None:
if is_primitive_axes_auto:
print("Primitive matrix (Auto):")
else:
print("Primitive matrix:")
for v in primitive_matrix:
print(f" {v}")
def print_supercell_matrix(
supercell_matrix: Union[Sequence, np.ndarray],
phonon_supercell_matrix: Optional[Union[Sequence, np.ndarray]] = None,
supercell_matrix: Sequence | np.ndarray,
phonon_supercell_matrix: Sequence | np.ndarray | None = None,
):
"""Print supercell matrix."""
if (np.diag(np.diag(supercell_matrix)) - supercell_matrix).any():

20
phono3py/cui/triplets_info.py
View File

@ -34,10 +34,12 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from collections.abc import Sequence
from typing import Optional, Union
from __future__ import annotations
import os
import numpy as np
from numpy.typing import ArrayLike
from phonopy.structure.atoms import PhonopyAtoms
from phono3py.file_IO import write_grid_address_to_hdf5, write_ir_grid_points
@ -48,13 +50,13 @@ from phono3py.phonon3.triplets import get_triplets_at_q
def write_grid_points(
primitive: PhonopyAtoms,
bz_grid: BZGrid,
band_indices: Optional[Union[Sequence, np.ndarray]] = None,
sigmas: Optional[Union[Sequence, np.ndarray]] = None,
temperatures: Optional[Union[Sequence, np.ndarray]] = None,
band_indices: ArrayLike | None = None,
sigmas: ArrayLike | None = None,
temperatures: ArrayLike | None = None,
is_kappa_star: bool = True,
is_lbte: bool = False,
compression: Union[str, int] = "gzip",
filename: bool = None,
compression: str | int = "gzip",
filename: str | os.PathLike | None = None,
):
"""Write grid points into files."""
ir_grid_points, ir_grid_weights = _get_ir_grid_points(
@ -109,8 +111,8 @@ def write_grid_points(
def show_num_triplets(
primitive: PhonopyAtoms,
bz_grid: BZGrid,
band_indices: Optional[Union[Sequence, np.ndarray]] = None,
grid_points: Optional[Union[Sequence, np.ndarray]] = None,
band_indices: ArrayLike | None = None,
grid_points: ArrayLike | None = None,
is_kappa_star: bool = True,
):
"""Show numbers of triplets at grid points."""

20
phono3py/file_IO.py
View File

@ -54,6 +54,7 @@ from phonopy.file_IO import (
write_FORCE_SETS,
)
from phono3py.phonon.grid import BZGrid
from phono3py.version import __version__
@ -483,9 +484,9 @@ def write_grid_address_to_hdf5(
grid_address,
mesh,
grid_mapping_table,
bz_grid=None,
compression: Union[str, int] = "gzip",
filename=None,
bz_grid: BZGrid | None = None,
compression: str | int = "gzip",
filename: str | os.PathLike | None = None,
):
"""Write grid addresses to grid_address.hdf5."""
suffix = _get_filename_suffix(mesh, filename=filename)
@ -1155,7 +1156,7 @@ def read_gamma_from_hdf5(
def read_collision_from_hdf5(
mesh,
indices=None,
indices: str | Sequence = "all",
grid_point=None,
band_index=None,
sigma=None,
@ -1643,7 +1644,11 @@ def parse_disp_fc3_yaml(filename="disp_fc3.yaml", return_cell=False):
return new_dataset
def parse_FORCES_FC2(disp_dataset, filename="FORCES_FC2", unit_conversion_factor=None):
def parse_FORCES_FC2(
disp_dataset: dict,
filename: str | os.PathLike = "FORCES_FC2",
unit_conversion_factor: float | None = None,
):
"""Parse type1 FORCES_FC2 file and store forces in disp_dataset."""
num_atom = disp_dataset["natom"]
num_disp = len(disp_dataset["first_atoms"])
@ -1664,7 +1669,10 @@ def parse_FORCES_FC2(disp_dataset, filename="FORCES_FC2", unit_conversion_factor
def parse_FORCES_FC3(
disp_dataset, filename="FORCES_FC3", use_loadtxt=False, unit_conversion_factor=None
disp_dataset: dict,
filename: str | os.PathLike = "FORCES_FC3",
use_loadtxt: bool = False,
unit_conversion_factor: float | None = None,
):
"""Parse type1 FORCES_FC3 and store forces in disp_dataset."""
num_atom = disp_dataset["natom"]

69
phono3py/interface/fc_calculator.py
View File

@ -116,14 +116,15 @@ class FC3Solver(FCSolver):
def _get_displacements_and_forces(self):
"""Return displacements and forces for fc3."""
assert self._dataset is not None
return get_displacements_and_forces_fc3(self._dataset)
def extract_fc2_fc3_calculators(
fc_calculator: Literal["traditional", "symfc", "alm"] | str | None,
order: int,
) -> Literal["traditional", "symfc", "alm"] | str | None:
"""Extract fc_calculator and fc_calculator_options for fc2 and fc3.
) -> Literal["traditional", "symfc", "alm"] | None:
"""Extract fc_calculator for fc2 and fc3.
fc_calculator : str
FC calculator. "|" separates fc2 and fc3. First and last
@ -133,20 +134,51 @@ def extract_fc2_fc3_calculators(
"""
if fc_calculator is None:
return fc_calculator
elif isinstance(fc_calculator, str):
if "|" in fc_calculator:
_fc_calculator = fc_calculator.split("|")[order - 2]
if _fc_calculator == "":
return None
else:
if fc_calculator.strip() == "":
return None
else:
_fc_calculator = fc_calculator
return _fc_calculator
return None
else:
raise RuntimeError("fc_calculator should be str or None.")
_fc_calculator = _split_fc_calculator_str(fc_calculator, order)
if _fc_calculator is None:
return None
fc_calculator_lower = _fc_calculator.lower()
if fc_calculator_lower not in ("traditional", "symfc", "alm"):
raise ValueError(
f"Unknown fc_calculator: {_fc_calculator}. "
"Available calculators are 'traditional', 'symfc', and 'alm'."
)
return fc_calculator_lower
def extract_fc2_fc3_calculators_options(
fc_calculator_opts: str | None,
order: int,
) -> str | None:
"""Extract fc_calculator_options for fc2 and fc3.
fc_calculator_opts : str
FC calculator options. "|" separates fc2 and fc3. First and last
parts separated correspond to fc2 and fc3 calculators, respectively.
order : int = 2 or 3
2 and 3 indicate fc2 and fc3, respectively.
"""
if fc_calculator_opts is None:
return None
else:
_fc_calculator_opts = _split_fc_calculator_str(fc_calculator_opts, order)
return _fc_calculator_opts
def _split_fc_calculator_str(fc_calculator: str, order: int) -> str | None:
if "|" in fc_calculator:
_fc_calculator = fc_calculator.split("|")[order - 2]
if _fc_calculator == "":
return None
else:
if fc_calculator.strip() == "":
return None
else:
_fc_calculator = fc_calculator
return _fc_calculator
def update_cutoff_fc_calculator_options(
@ -235,6 +267,11 @@ def determine_cutoff_pair_distance(
"available for symfc calculator."
)
symfc_options = {"memsize": {3: _symfc_memory_size}}
if supercell is None or primitive is None or symmetry is None:
raise RuntimeError(
"supercell, primitive, and symmetry are required to estimate "
"cutoff_pair_distance by memory size."
)
update_symfc_cutoff_by_memsize(
symfc_options, supercell, primitive, symmetry, verbose=log_level > 0
)
@ -283,7 +320,7 @@ def _get_cutoff_pair_distance(
cutoff_pair_distance,
)
symfc_options = parse_symfc_options(
extract_fc2_fc3_calculators(_fc_calculator_options, 3), 3
extract_fc2_fc3_calculators_options(_fc_calculator_options, 3), 3
)
_cutoff_pair_distance = cutoff_pair_distance

25
phono3py/other/isotope.py
View File

@ -36,9 +36,8 @@
from __future__ import annotations
from typing import Optional, Union
import numpy as np
from numpy.typing import ArrayLike
from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
from phonopy.phonon.tetrahedron_mesh import get_tetrahedra_frequencies
from phonopy.physical_units import get_physical_units
@ -58,9 +57,9 @@ from phono3py.phonon.solver import run_phonon_solver_c, run_phonon_solver_py
def get_mass_variances(
primitive: Optional[PhonopyAtoms] = None,
symbols: Optional[Union[list[str], tuple[str]]] = None,
isotope_data: Optional[dict] = None,
primitive: PhonopyAtoms | None = None,
symbols: list[str] | tuple[str] | None = None,
isotope_data: dict | None = None,
):
"""Calculate mass variances."""
if primitive is not None:
@ -93,14 +92,14 @@ class Isotope:
def __init__(
self,
mesh,
primitive,
mesh: float | ArrayLike,
primitive: Primitive,
mass_variances=None, # length of list is num_atom.
isotope_data=None,
band_indices=None,
sigma=None,
bz_grid=None,
frequency_factor_to_THz=None,
bz_grid: BZGrid | None = None,
frequency_factor_to_THz: float | None = None,
use_grg=False,
symprec=1e-5,
cutoff_frequency=None,
@ -116,7 +115,6 @@ class Isotope:
self._mass_variances = np.array(mass_variances, dtype="double")
self._primitive = primitive
self._sigma = sigma
self._bz_grid = bz_grid
self._symprec = symprec
if cutoff_frequency is None:
self._cutoff_frequency = 0
@ -143,7 +141,7 @@ class Isotope:
else:
self._band_indices = np.array(band_indices, dtype="int64")
if self._bz_grid is None:
if bz_grid is None:
primitive_symmetry = Symmetry(self._primitive, self._symprec)
self._bz_grid = BZGrid(
self._mesh,
@ -152,6 +150,8 @@ class Isotope:
use_grg=use_grg,
store_dense_gp_map=True,
)
else:
self._bz_grid = bz_grid
def set_grid_point(self, grid_point):
"""Initialize grid points."""
@ -196,8 +196,9 @@ class Isotope:
return self._gamma
@property
def bz_grid(self):
def bz_grid(self) -> BZGrid:
"""Return BZgrid class instance."""
assert self._bz_grid is not None
return self._bz_grid
@property

12
phono3py/phonon/grid.py
View File

@ -144,12 +144,12 @@ class BZGrid:
def __init__(
self,
mesh: Union[int, float, Sequence, np.ndarray],
reciprocal_lattice=None,
lattice=None,
symmetry_dataset: Optional[SpglibDataset] = None,
transformation_matrix: Optional[Union[Sequence, np.ndarray]] = None,
is_shift: Optional[Union[list, np.ndarray]] = None,
mesh: float | ArrayLike,
reciprocal_lattice: ArrayLike | None = None,
lattice: ArrayLike | None = None,
symmetry_dataset: SpglibDataset | None = None,
transformation_matrix: ArrayLike | None = None,
is_shift: ArrayLike | None = None,
is_time_reversal: bool = True,
use_grg: bool = False,
force_SNF: bool = False,

24
phono3py/phonon/solver.py
View File

@ -34,7 +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 numpy as np
from phonopy.harmonic.dynamical_matrix import (
DynamicalMatrix,
DynamicalMatrixGL,
DynamicalMatrixNAC,
)
from phonopy.physical_units import get_physical_units
from phonopy.structure.cells import sparse_to_dense_svecs
@ -81,8 +88,8 @@ def run_phonon_solver_c(
'U' or 'L' for lapack zheev solver. Default is 'L'.
"""
import phono3py._phono3py as phono3c
import phono3py._phononcalc as phononcalc
import phono3py._phono3py as phono3c # type: ignore[import-untyped]
import phono3py._phononcalc as phononcalc # type: ignore[import-untyped]
if frequency_conversion_factor is None:
_frequency_conversion_factor = get_physical_units().DefaultToTHz
@ -100,7 +107,7 @@ def run_phonon_solver_c(
dielectric,
) = _extract_params(dm)
if dm.is_nac() and dm.nac_method == "gonze":
if isinstance(dm, DynamicalMatrixGL):
gonze_nac_dataset = dm.Gonze_nac_dataset
if gonze_nac_dataset[0] is None:
dm.make_Gonze_nac_dataset()
@ -112,6 +119,7 @@ def run_phonon_solver_c(
G_list, # List of G points where d-d interactions are integrated.
Lambda,
) = gonze_nac_dataset # Convergence parameter
assert Lambda is not None
fc = gonze_fc
use_GL_NAC = True
else:
@ -120,7 +128,7 @@ def run_phonon_solver_c(
dd_q0 = np.zeros(2) # dummy variable
G_list = np.zeros(3) # dummy variable
Lambda = 0 # dummy variable
if not dm.is_nac():
if not isinstance(dm, DynamicalMatrixNAC):
born = np.zeros((3, 3)) # dummy variable
dielectric = np.zeros(3) # dummy variable
fc = dm.force_constants
@ -168,7 +176,7 @@ def run_phonon_solver_c(
dd_q0,
G_list,
float(Lambda),
dm.is_nac() * 1,
isinstance(dm, DynamicalMatrixNAC) * 1,
is_nac_q_zero * 1,
use_GL_NAC * 1,
lapack_zheev_uplo,
@ -207,14 +215,14 @@ def run_phonon_solver_py(
dynamical_matrix.run(q)
dm = dynamical_matrix.dynamical_matrix
eigvals, eigvecs = np.linalg.eigh(dm, UPLO=lapack_zheev_uplo)
eigvals = eigvals.real
eigvals = eigvals.real # type: ignore[no-untyped-call]
frequencies[gp] = (
np.sqrt(np.abs(eigvals)) * np.sign(eigvals) * frequency_conversion_factor
)
eigenvectors[gp] = eigvecs
def _extract_params(dm):
def _extract_params(dm: DynamicalMatrix | DynamicalMatrixNAC):
svecs, multi = dm.primitive.get_smallest_vectors()
if dm.primitive.store_dense_svecs:
_svecs = svecs
@ -225,7 +233,7 @@ def _extract_params(dm):
masses = np.array(dm.primitive.masses, dtype="double")
rec_lattice = np.array(np.linalg.inv(dm.primitive.cell), dtype="double", order="C")
positions = np.array(dm.primitive.positions, dtype="double", order="C")
if dm.is_nac():
if isinstance(dm, DynamicalMatrixNAC):
born = dm.born
nac_factor = dm.nac_factor
dielectric = dm.dielectric_constant

7
phono3py/phonon/velocity_operator.py
View File

@ -35,6 +35,7 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.harmonic.dynamical_matrix import DynamicalMatrixGL
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.physical_units import get_physical_units
@ -210,11 +211,7 @@ class VelocityOperator(GroupVelocity):
if np.linalg.norm(q) < np.linalg.norm(delta_q):
flag_gamma = True
if (
(self._dynmat.is_nac())
and (self._dynmat.nac_method == "gonze")
and flag_gamma
):
if isinstance(dynmat, DynamicalMatrixGL) and flag_gamma:
dynmat.run(
q - delta_q, q_direction=(q - delta_q) / np.linalg.norm(q - delta_q)
)

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.

57
phono3py/phonon3/fc3.py
View File

@ -34,10 +34,13 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import logging
import sys
import numpy as np
from numpy.typing import NDArray
from phonopy.harmonic.force_constants import (
distribute_force_constants,
get_nsym_list_and_s2pp,
@ -573,30 +576,50 @@ def cutoff_fc3_by_zero(fc3, supercell, cutoff_distance, p2s_map=None, symprec=1e
break
def show_drift_fc3(fc3, primitive=None, name="fc3"):
"""Show drift of fc3."""
def show_drift_fc3(
fc3: NDArray, primitive: Primitive | None = None, name: str = "fc3", digit: int = 8
):
"""Show max drift of fc3."""
maxval1, maxval2, maxval3, xyz1, xyz2, xyz3 = get_drift_fc3(
fc3, primitive=primitive
)
text = f"Max drift of {name}: "
text += f"{maxval1:.{digit}f} ({'xyz'[xyz1[0]]}{'xyz'[xyz1[1]]}{'xyz'[xyz1[2]]}) "
text += f"{maxval2:.{digit}f} ({'xyz'[xyz2[0]]}{'xyz'[xyz2[1]]}{'xyz'[xyz2[2]]}) "
text += f"{maxval3:.{digit}f} ({'xyz'[xyz3[0]]}{'xyz'[xyz3[1]]}{'xyz'[xyz3[2]]})"
print(text)
def get_drift_fc3(
fc3: NDArray, primitive: Primitive | None = None
) -> tuple[float, float, float, list[int], list[int], list[int]]:
"""Return max drift of fc3."""
if fc3.shape[0] == fc3.shape[1]:
num_atom = fc3.shape[0]
maxval1 = 0
maxval2 = 0
maxval3 = 0
klm1 = [0, 0, 0]
klm2 = [0, 0, 0]
klm3 = [0, 0, 0]
xyz1 = [0, 0, 0]
xyz2 = [0, 0, 0]
xyz3 = [0, 0, 0]
for i, j, k, ll, m in list(np.ndindex((num_atom, num_atom, 3, 3, 3))):
val1 = fc3[:, i, j, k, ll, m].sum()
val2 = fc3[i, :, j, k, ll, m].sum()
val3 = fc3[i, j, :, k, ll, m].sum()
if abs(val1) > abs(maxval1):
maxval1 = val1
klm1 = [k, ll, m]
xyz1 = [k, ll, m]
if abs(val2) > abs(maxval2):
maxval2 = val2
klm2 = [k, ll, m]
xyz2 = [k, ll, m]
if abs(val3) > abs(maxval3):
maxval3 = val3
klm3 = [k, ll, m]
xyz3 = [k, ll, m]
else:
if primitive is None:
raise RuntimeError(
"Primitive cell is required to get drift of compact fc3."
)
try:
import phono3py._phono3py as phono3c
@ -614,9 +637,9 @@ def show_drift_fc3(fc3, primitive=None, name="fc3"):
maxval1 = 0
maxval2 = 0
maxval3 = 0
klm1 = [0, 0, 0]
klm2 = [0, 0, 0]
klm3 = [0, 0, 0]
xyz1 = [0, 0, 0]
xyz2 = [0, 0, 0]
xyz3 = [0, 0, 0]
phono3c.transpose_compact_fc3(
fc3, permutations, s2pp_map, p2s_map, nsym_list, 0
) # dim[0] <--> dim[1]
@ -624,7 +647,7 @@ def show_drift_fc3(fc3, primitive=None, name="fc3"):
val1 = fc3[i, :, j, k, ll, m].sum()
if abs(val1) > abs(maxval1):
maxval1 = val1
klm1 = [k, ll, m]
xyz1 = [k, ll, m]
phono3c.transpose_compact_fc3(
fc3, permutations, s2pp_map, p2s_map, nsym_list, 0
) # dim[0] <--> dim[1]
@ -633,10 +656,10 @@ def show_drift_fc3(fc3, primitive=None, name="fc3"):
val3 = fc3[i, j, :, k, ll, m].sum()
if abs(val2) > abs(maxval2):
maxval2 = val2
klm2 = [k, ll, m]
xyz2 = [k, ll, m]
if abs(val3) > abs(maxval3):
maxval3 = val3
klm3 = [k, ll, m]
xyz3 = [k, ll, m]
except ImportError as exc:
text = (
"Import error at phono3c.tranpose_compact_fc3. "
@ -644,11 +667,7 @@ def show_drift_fc3(fc3, primitive=None, name="fc3"):
)
raise RuntimeError(text) from exc
text = "Max drift of %s: " % name
text += "%f (%s%s%s) " % (maxval1, "xyz"[klm1[0]], "xyz"[klm1[1]], "xyz"[klm1[2]])
text += "%f (%s%s%s) " % (maxval2, "xyz"[klm2[0]], "xyz"[klm2[1]], "xyz"[klm2[2]])
text += "%f (%s%s%s)" % (maxval3, "xyz"[klm3[0]], "xyz"[klm3[1]], "xyz"[klm3[2]])
print(text)
return maxval1, maxval2, maxval3, xyz1, xyz2, xyz3
def _set_permutation_symmetry_fc3_elem_with_cutoff(fc3, fc3_done, a, b, c):

12
phono3py/phonon3/gruneisen.py
View File

@ -37,7 +37,11 @@
import sys
import numpy as np
from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
from phonopy.harmonic.dynamical_matrix import (
DynamicalMatrixGL,
DynamicalMatrixNAC,
get_dynamical_matrix,
)
from phonopy.physical_units import get_physical_units
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import Primitive
@ -96,8 +100,8 @@ def run_gruneisen_parameters(
if log_level > 0:
dm = gruneisen.dynamical_matrix
if dm.is_nac() and dm.nac_method == "gonze":
dm.show_Gonze_nac_message()
if isinstance(dm, DynamicalMatrixGL):
dm.show_nac_message()
if mesh is not None:
gruneisen.set_sampling_mesh(mesh, rotations=rotations, is_gamma_center=True)
@ -310,7 +314,7 @@ class Gruneisen:
gruneisen_parameters = []
frequencies = []
for i, q in enumerate(qpoints):
if self._dm.is_nac():
if isinstance(self._dm, DynamicalMatrixNAC):
if (np.abs(q) < 1e-5).all(): # If q is almost at Gamma
if self._run_mode == "band":
# Direction estimated from neighboring point

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(

2
phono3py/version.py
View File

@ -34,4 +34,4 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
__version__ = "3.17.1"
__version__ = "3.18.0"

2
pyproject.toml
View File

@ -16,7 +16,7 @@ dependencies = [
"matplotlib",
"h5py",
"spglib",
"phonopy>=2.41,<2.42",
"phonopy>=2.42,<2.43",
]
license = "BSD-3-Clause"
license-files = ["LICENSE"]

62
test/api/test_api_phono3py.py
View File

@ -6,9 +6,11 @@ from pathlib import Path
import numpy as np
import pytest
from phonopy.harmonic.force_constants import get_drift_force_constants
from phonopy.interface.pypolymlp import PypolymlpParams
from phono3py import Phono3py
from phono3py.phonon3.fc3 import get_drift_fc3
cwd = Path(__file__).parent
@ -207,3 +209,63 @@ def test_use_pypolymlp_mgo(mgo_222rd_444rd: Phono3py):
assert (
pytest.approx(63.0018546, abs=1e-1) == ph3.thermal_conductivity.kappa[0, 0, 0]
)
@pytest.mark.parametrize("is_compact_fc", [True, False])
def test_symmetrize_fc_traditional(si_pbesol: Phono3py, is_compact_fc: bool):
"""Test symmetrize_fc3 and symmetrize_fc2 with traditional approach."""
ph3 = Phono3py(
si_pbesol.unitcell,
supercell_matrix=si_pbesol.supercell_matrix,
primitive_matrix=si_pbesol.primitive_matrix,
log_level=2,
)
ph3.dataset = si_pbesol.dataset
ph3.produce_fc3(is_compact_fc=is_compact_fc)
assert ph3.fc3 is not None
assert ph3.fc2 is not None
v1, v2, v3, _, _, _ = get_drift_fc3(ph3.fc3, primitive=ph3.primitive)
np.testing.assert_allclose(
np.abs([v1, v2, v3]), [1.755065e-01, 1.749287e-01, 3.333333e-05], atol=1e-6
)
ph3.symmetrize_fc3(options="level=3")
v1_sym, v2_sym, v3_sym, _, _, _ = get_drift_fc3(ph3.fc3, primitive=ph3.primitive)
if is_compact_fc:
np.testing.assert_allclose(
np.abs([v1_sym, v2_sym, v3_sym]), 1.217081e-05, atol=1e-6
)
else:
np.testing.assert_allclose(
np.abs([v1_sym, v2_sym, v3_sym]), 1.421085e-14, atol=1e-6
)
v1_sym, v2_sym, _, _ = get_drift_force_constants(ph3.fc2, primitive=ph3.primitive)
if is_compact_fc:
np.testing.assert_allclose(np.abs([v1_sym, v2_sym]), 1.0e-06, atol=1e-6)
else:
np.testing.assert_allclose(np.abs([v1_sym, v2_sym]), 1.0e-06, atol=1e-6)
@pytest.mark.parametrize("is_compact_fc", [True, False])
def test_symmetrize_fc_symfc(si_pbesol: Phono3py, is_compact_fc: bool):
"""Test symmetrize_fc3 and symmetrize_fc2 with symfc."""
pytest.importorskip("symfc")
ph3 = Phono3py(
si_pbesol.unitcell,
supercell_matrix=si_pbesol.supercell_matrix,
primitive_matrix=si_pbesol.primitive_matrix,
log_level=2,
)
ph3.dataset = si_pbesol.dataset
ph3.produce_fc3(is_compact_fc=is_compact_fc)
assert ph3.fc3 is not None
assert ph3.fc2 is not None
ph3.symmetrize_fc3(use_symfc_projector=True)
v1_sym, v2_sym, v3_sym, _, _, _ = get_drift_fc3(ph3.fc3, primitive=ph3.primitive)
np.testing.assert_allclose([v1_sym, v2_sym, v3_sym], 0, atol=1e-6)
ph3.symmetrize_fc2(use_symfc_projector=True)
v1_sym, v2_sym, _, _ = get_drift_force_constants(ph3.fc2, primitive=ph3.primitive)
np.testing.assert_allclose([v1_sym, v2_sym], 0, atol=1e-6)

2
test/phonon3/test_fc3.py
View File

@ -72,7 +72,7 @@ def test_cutoff_fc3_zero_compact_fc(nacl_pbe_compact_fc: Phono3py):
fc3 = ph.fc3.copy()
cutoff_fc3_by_zero(fc3, ph.supercell, 5, p2s_map=ph.primitive.p2s_map)
abs_delta = np.abs(ph.fc3 - fc3).sum()
assert np.abs(164.359250 - abs_delta) < 1e-3
assert np.abs(164.350663 - abs_delta) < 1e-3
def test_fc3(si_pbesol_111: Phono3py):

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)