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

Implement relaxation using pypolymlp in CUI

This commit is contained in:
Atsushi Togo 2025-07-20 18:42:07 +09:00
parent cbcf9e84e5
commit 3564186adc
10 changed files with 246 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
doc/changelog.md
View File

@ -5,6 +5,8 @@
## Jul-xx-2025: Version 3.17.2
- 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

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.

85
phono3py/api_phono3py.py
View File

@ -847,7 +847,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
@ -1558,10 +1558,10 @@ class Phono3py:
if symmetrize_fc3r:
if use_symfc_projector and fc_calculator is None:
self.symmetrize_fc3(
use_symfc_projector=True, symfc_options=fc_calculator_options
use_symfc_projector=True, options=fc_calculator_options
)
else:
self.symmetrize_fc3()
self.symmetrize_fc3(options=fc_calculator_options)
elif fc_calculator == "symfc":
symfc_solver = cast(SymfcFCSolver, fc_solver.fc_solver)
fc3_nonzero_elems = symfc_solver.get_nonzero_atomic_indices_fc3()
@ -1587,30 +1587,29 @@ class Phono3py:
if use_symfc_projector and fc_calculator is None:
self.symmetrize_fc2(
use_symfc_projector=True,
symfc_options=fc_calculator_options,
options=fc_calculator_options,
)
else:
self.symmetrize_fc2()
self.symmetrize_fc2(options=fc_calculator_options)
def symmetrize_fc3(
self,
level: int = 1,
use_symfc_projector: bool = False,
symfc_options: str | None = None,
options: str | None = None,
):
"""Symmetrize fc3 by symfc projector or traditional approach.
Parameters
----------
level : int, optional
Number of times translational and permutation symmetries are applied
for traditional symmetrization. Default is 1.
use_symfc_projector : bool, optional
If True, the force constants are symmetrized by symfc projector
instead of traditional approach.
symfc_options : str or None, optional
Options for symfc projector. "use_mkl=true" calls sparse_dot_mkl
(required to install it).
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:
@ -1619,22 +1618,33 @@ class Phono3py:
if use_symfc_projector:
if self._log_level:
print("Symmetrizing fc3 by symfc projector.", flush=True)
if symfc_options is None:
options = None
if options is None:
_options = None
else:
options = parse_symfc_options(symfc_options, 3)
_options = parse_symfc_options(options, 3)
self._fc3 = symmetrize_by_projector(
self._supercell,
self._fc3,
3,
primitive=self._primitive,
options=options,
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)
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)
@ -1703,28 +1713,28 @@ 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,
level: int = 1,
use_symfc_projector: bool = False,
symfc_options: str | None = None,
options: str | None = None,
):
"""Symmetrize fc2 by symfc projector or traditional approach.
Parameters
----------
level : int, optional
Number of times translational and permutation symmetries are applied
for traditional symmetrization. Default is 1.
use_symfc_projector : bool, optional
If True, the force constants are symmetrized by symfc projector
instead of traditional approach.
symfc_options : str or None, optional
Options for symfc projector. "use_mkl=true" calls sparse_dot_mkl
(required to install it).
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:
@ -1743,22 +1753,33 @@ class Phono3py:
if use_symfc_projector:
if self._log_level:
print("Symmetrizing fc2 by symfc projector.", flush=True)
if symfc_options is None:
options = None
if options is None:
_options = None
else:
options = parse_symfc_options(symfc_options, 2)
_options = parse_symfc_options(options, 2)
self._fc2 = symmetrize_by_projector(
supercell,
self._fc2,
2,
primitive=primitive,
options=options,
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)
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)

6
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:

7
phono3py/cui/phono3py_argparse.py
View File

@ -727,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",

130
phono3py/cui/phono3py_script.py
View File

@ -57,14 +57,20 @@ 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
@ -431,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,
@ -439,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,
@ -475,9 +507,7 @@ 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(
@ -563,7 +593,7 @@ def _run_pypolymlp(
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
@ -575,9 +605,57 @@ def _run_pypolymlp(
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 as e:
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,
_ph3py,
displacement_distance=settings.displacement_distance,
number_of_snapshots=settings.random_displacements,
number_estimation_factor=settings.rd_number_estimation_factor,
@ -590,7 +668,7 @@ def _run_pypolymlp(
)
if log_level:
print(f'Dataset generated using MLPs was written in "{mlp_eval_filename}".')
ph3py.save(mlp_eval_filename)
_ph3py.save(mlp_eval_filename)
else:
if log_level:
print(
@ -598,9 +676,11 @@ def _run_pypolymlp(
"calculations."
)
_finalize_phono3py(
ph3py, confs_dict, log_level=log_level, filename=output_yaml_filename
_ph3py, confs_dict, log_level=log_level, filename=output_yaml_filename
)
return _ph3py
def _produce_force_constants(
ph3py: Phono3py,
@ -1096,7 +1176,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
)
@ -1183,6 +1263,7 @@ def main(**argparse_control):
####################
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
@ -1266,7 +1347,8 @@ def main(**argparse_control):
if ph3py.fc3 is None or (
ph3py.fc2 is None and ph3py.phonon_supercell_matrix is None
):
_run_pypolymlp(
# Note that ph3py is replaced when relax_atomic_positions is True.
ph3py = _run_pypolymlp(
ph3py,
settings,
confs_dict,

11
phono3py/file_IO.py
View File

@ -1644,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"])
@ -1665,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"]

2
test/api/test_api_phono3py.py
View File

@ -229,7 +229,7 @@ def test_symmetrize_fc_traditional(si_pbesol: Phono3py, is_compact_fc: bool):
np.testing.assert_allclose(
np.abs([v1, v2, v3]), [1.755065e-01, 1.749287e-01, 3.333333e-05], atol=1e-6
)
ph3.symmetrize_fc3(level=3)
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(

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):