Update pypolymlp feature

phono3py/api_phono3py.py

@@ -40,7 +40,6 @@ from collections.abc import Sequence
 from typing import Literal, Optional, Union
 
 import numpy as np
-from phonopy import Phonopy
 from phonopy.exception import ForceCalculatorRequiredError
 from phonopy.harmonic.displacement import (
     directions_to_displacement_dataset,
@@ -56,7 +55,13 @@ from phonopy.harmonic.force_constants import (
     symmetrize_force_constants,
 )
 from phonopy.interface.fc_calculator import get_fc2
-from phonopy.interface.pypolymlp import PypolymlpParams
+from phonopy.interface.pypolymlp import (
+    PypolymlpData,
+    PypolymlpParams,
+    develop_polymlp,
+    evalulate_polymlp,
+    parse_mlp_params,
+)
 from phonopy.structure.atoms import PhonopyAtoms
 from phonopy.structure.cells import (
     Primitive,
@@ -295,8 +300,9 @@ class Phono3py:
         self._fc2 = None
         self._fc3 = None
 
-        # MLP container
+        # MLP
-        self._mlp_phonopy = None
+        self._mlp = None
+        self._mlp_dataset = None
 
         # Setup interaction
         self._interaction = None
@@ -681,6 +687,35 @@ class Phono3py:
 
         self._phonon_supercells_with_displacements = None
 
+    @property
+    def mlp_dataset(self) -> Optional[dict]:
+        """Return displacement-force dataset.
+
+        The supercell matrix is equal to that of usual displacement-force
+        dataset. Only type 2 format is supported. "displacements",
+        "forces", and "supercell_energies" should be contained.
+
+        """
+        return self._mlp_dataset
+
+    @mlp_dataset.setter
+    def mlp_dataset(self, mlp_dataset: dict):
+        if not isinstance(mlp_dataset, dict):
+            raise TypeError("mlp_dataset has to be a dictionary.")
+        if "displacements" not in mlp_dataset:
+            raise RuntimeError("Displacements have to be given.")
+        if "forces" not in mlp_dataset:
+            raise RuntimeError("Forces have to be given.")
+        if "supercell_energy" in mlp_dataset:
+            raise RuntimeError("Supercell energies have to be given.")
+        if len(mlp_dataset["displacements"]) != len(mlp_dataset["forces"]):
+            raise RuntimeError("Length of displacements and forces are different.")
+        if len(mlp_dataset["displacements"]) != len(mlp_dataset["supercell_energies"]):
+            raise RuntimeError(
+                "Length of displacements and supercell_energies are different."
+            )
+        self._mlp_dataset = mlp_dataset
+
     @property
     def band_indices(self) -> list[np.ndarray]:
         """Setter and getter of band indices.
@@ -811,7 +846,7 @@ class Phono3py:
             for disp1 in dataset["first_atoms"]:
                 num_scells += len(disp1["second_atoms"])
             displacements = np.zeros(
-                (num_scells, self._supercell.get_number_of_atoms(), 3),
+                (num_scells, len(self._supercell), 3),
                 dtype="double",
                 order="C",
             )
@@ -990,30 +1025,6 @@ class Phono3py:
     def phonon_supercell_energies(self, values):
         self._set_phonon_forces_energies(values, target="supercell_energies")
 
-    @property
-    def mlp_dataset(self) -> Optional[dict]:
-        """Return displacement-force dataset.
-
-        The supercell matrix is equal to that of usual displacement-force
-        dataset. Only type 2 format is supported. "displacements",
-        "forces", and "supercell_energies" should be contained.
-
-        """
-        if self._mlp_phonopy is None:
-            return None
-        return self._mlp_phonopy.mlp_dataset
-
-    @mlp_dataset.setter
-    def mlp_dataset(self, mlp_dataset: dict):
-        if self._mlp_phonopy is None:
-            self._mlp_phonopy = Phonopy(
-                self._unitcell,
-                supercell_matrix=self._supercell_matrix,
-                symprec=self._symprec,
-                log_level=self._log_level,
-            )
-            self._mlp_phonopy.mlp_dataset = mlp_dataset
-
     @property
     def phph_interaction(self):
         """Return Interaction instance."""
@@ -1385,53 +1396,31 @@ class Phono3py:
         symmetrize_fc3r: bool = False,
         is_compact_fc: bool = False,
         fc_calculator: Optional[str] = None,
-        fc_calculator_options: Optional[str] = None,
+        fc_calculator_options: Optional[Union[str, dict]] = None,
-        use_pypolymlp: bool = False,
-        mlp_params: Optional[Union[PypolymlpParams, dict]] = None,
     ):
         """Calculate fc3 from displacements and forces.
 
         Parameters
         ----------
-        symmetrize_fc3r : bool
+        symmetrize_fc3r : bool, optional
             Only for type 1 displacement_dataset, translational and permutation
             symmetries are applied after creating fc3. This symmetrization is
             not very sophisticated and can break space group symmetry, but often
             useful. If better symmetrization is expected, it is recommended to
             use external force constants calculator such as ALM. Default is
             False.
-        is_compact_fc : bool
+        is_compact_fc : bool, optional
             fc3 shape is
                 False: (supercell, supercell, supecell, 3, 3, 3) True:
                 (primitive, supercell, supecell, 3, 3, 3)
             where 'supercell' and 'primitive' indicate number of atoms in these
             cells. Default is False.
-        fc_calculator : str or None
+        fc_calculator : str, optional
             Force constants calculator given by str.
-        fc_calculator_options : dict
+        fc_calculator_options : dict or str, optional
             Options for external force constants calculator.
-        use_pypolymlp : bool
-            Use MLP of pypolymlp to calculate fc3. Default is False. mlp_dataset
-            has to be set before calling this method.
-        mlp_params : PypolymlpParams or dict, optional
-            Parameters for developing MLP. Default is None. When dict is given,
-            PypolymlpParams instance is created from the dict.
 
         """
-        if use_pypolymlp:
-            if self._mlp_phonopy is None:
-                msg = "mlp_dataset has to be set before calling this method."
-                raise RuntimeError(msg)
-            if self._dataset is None or "displacements" not in self._dataset:
-                raise RuntimeError(
-                    "Type 2 displacements are not set. Run generate_displacements."
-                )
-            self._mlp_phonopy.develop_mlp(params=mlp_params)
-            self._mlp_phonopy.displacements = self._dataset["displacements"]
-            self._mlp_phonopy.evaluate_mlp()
-            self._dataset["forces"] = self._mlp_phonopy.forces
-            self._dataset["supercell_energies"] = self._mlp_phonopy.supercell_energies
-
         fc3_calculator, fc3_calculator_options = self._extract_fc2_fc3_calculators(
             fc_calculator, fc_calculator_options, 3
         )
@@ -2172,6 +2161,65 @@ class Phono3py:
         with open(filename, "w") as w:
             w.write(str(ph3py_yaml))
 
+    def develop_mlp(self, params: Optional[Union[PypolymlpParams, dict, str]] = None):
+        """Develop MLP of pypolymlp.
+
+        Parameters
+        ----------
+        params : PypolymlpParams or dict, optional
+            Parameters for developing MLP. Default is None. When dict is given,
+            PypolymlpParams instance is created from the dict.
+
+        """
+        if self._mlp_dataset is None:
+            raise RuntimeError("MLP dataset is not set.")
+
+        if params is not None:
+            _params = parse_mlp_params(params)
+
+        disps = self._mlp_dataset["displacements"]
+        forces = self._mlp_dataset["forces"]
+        energies = self._mlp_dataset["supercell_energies"]
+        n = int(len(disps) * 0.9)
+        train_data = PypolymlpData(
+            displacements=disps[:n], forces=forces[:n], supercell_energies=energies[:n]
+        )
+        test_data = PypolymlpData(
+            displacements=disps[n:], forces=forces[n:], supercell_energies=energies[n:]
+        )
+        self._mlp = develop_polymlp(
+            self._supercell,
+            train_data,
+            test_data,
+            params=_params,
+            verbose=self._log_level - 1 > 0,
+        )
+
+    def evaluate_mlp(self):
+        """Evaluate the machine learning potential of pypolymlp.
+
+        This method calculates the supercell energies and forces from the MLP
+        for the displacements in self._dataset of type 2. The results are stored
+        in self._dataset.
+
+        The displacements may be generated by the produce_force_constants method
+        with number_of_snapshots > 0. With MLP, a small distance parameter, such
+        as 0.001, can be numerically stable for the computation of force
+        constants.
+
+        """
+        if self._mlp is None:
+            raise RuntimeError("MLP is not developed yet.")
+
+        if self.supercells_with_displacements is None:
+            raise RuntimeError("Displacements are not set. Run generate_displacements.")
+
+        energies, forces, _ = evalulate_polymlp(
+            self._mlp, self.supercells_with_displacements
+        )
+        self.supercell_energies = energies
+        self.forces = forces
+
     ###################
     # private methods #
     ###################
@@ -2252,15 +2300,28 @@ class Phono3py:
                 raise RuntimeError
 
     def _build_phonon_supercells_with_displacements(
-        self, supercell: PhonopyAtoms, displacement_dataset
+        self, supercell: PhonopyAtoms, dataset
     ):
         supercells = []
+        positions = supercell.positions
         magmoms = supercell.magnetic_moments
         masses = supercell.masses
         numbers = supercell.numbers
         lattice = supercell.cell
 
-        for disp1 in displacement_dataset["first_atoms"]:
+        if "displacements" in dataset:
+            for disp in dataset["displacements"]:
+                supercells.append(
+                    PhonopyAtoms(
+                        numbers=numbers,
+                        masses=masses,
+                        magnetic_moments=magmoms,
+                        positions=positions + disp,
+                        cell=lattice,
+                    )
+                )
+        else:
+            for disp1 in dataset["first_atoms"]:
                 disp_cart1 = disp1["displacement"]
                 positions = supercell.positions
                 positions[disp1["number"]] += disp_cart1
@@ -2277,7 +2338,6 @@ class Phono3py:
         return supercells
 
     def _build_supercells_with_displacements(self):
-        supercells = []
         magmoms = self._supercell.magnetic_moments
         masses = self._supercell.masses
         numbers = self._supercell.numbers
@@ -2287,6 +2347,7 @@ class Phono3py:
             self._supercell, self._dataset
         )
 
+        if "first_atoms" in self._dataset:
             for disp1 in self._dataset["first_atoms"]:
                 disp_cart1 = disp1["displacement"]
                 for disp2 in disp1["second_atoms"]:

phono3py/cui/create_force_constants.py

@@ -40,7 +40,8 @@ import copy
 import os
 import pathlib
 import sys
-from typing import Optional
+from dataclasses import asdict
+from typing import Optional, Union
 
 import numpy as np
 from phonopy.cui.phonopy_script import file_exists, print_error
@@ -52,6 +53,7 @@ from phonopy.harmonic.force_constants import (
 )
 from phonopy.interface.calculator import get_default_physical_units
 from phonopy.interface.fc_calculator import fc_calculator_names
+from phonopy.interface.pypolymlp import PypolymlpParams, parse_mlp_params
 
 from phono3py import Phono3py
 from phono3py.cui.show_log import show_phono3py_force_constants_settings
@@ -122,6 +124,11 @@ def create_phono3py_force_constants(
                 settings.cutoff_pair_distance,
                 fc_calculator,
                 fc_calculator_options,
+                settings.use_pypolymlp,
+                settings.mlp_params,
+                settings.displacement_distance,
+                settings.random_displacements,
+                settings.random_seed,
                 log_level,
             )
 
@@ -249,7 +256,7 @@ def parse_forces(
     # Try to read FORCES_FC* if type-2 and return dataset.
     # None is returned unless type-2.
     # can emit FileNotFoundError.
-    if dataset is None or dataset is not None and not forces_in_dataset(dataset):
+    if dataset is None or (dataset is not None and not forces_in_dataset(dataset)):
         _dataset = read_type2_dataset(
             natom, filename=force_filename, log_level=log_level
         )
@@ -321,8 +328,10 @@ def forces_in_dataset(dataset: dict) -> bool:
     )
 
 
-def displacements_in_dataset(dataset: dict) -> bool:
+def displacements_in_dataset(dataset: Optional[dict]) -> bool:
     """Return whether displacements in dataset or not."""
+    if dataset is None:
+        return False
     return "displacements" in dataset or "first_atoms" in dataset
 
 
@@ -443,6 +452,11 @@ def _create_phono3py_fc3(
     cutoff_pair_distance: Optional[float],
     fc_calculator: Optional[str],
     fc_calculator_options: Optional[str],
+    use_pypoymplp: bool,
+    mlp_params: Union[str, dict, PypolymlpParams],
+    displacement_distance: Optional[float],
+    number_of_snapshots: Optional[int],
+    random_seed: Optional[int],
     log_level: int,
 ):
     """Read or calculate fc3.
@@ -481,6 +495,17 @@ def _create_phono3py_fc3(
         # from _get_type2_dataset
         file_exists(e.filename, log_level)
 
+    if use_pypoymplp:
+        phono3py.mlp_dataset = dataset
+        run_pypolymlp_to_compute_forces(
+            phono3py,
+            mlp_params=mlp_params,
+            displacement_distance=displacement_distance,
+            number_of_snapshots=number_of_snapshots,
+            random_seed=random_seed,
+            log_level=log_level,
+        )
+    else:
         phono3py.dataset = dataset
     phono3py.produce_fc3(
         symmetrize_fc3r=symmetrize_fc3r,
@@ -490,6 +515,68 @@ def _create_phono3py_fc3(
     )
 
 
+def run_pypolymlp_to_compute_forces(
+    ph3py: Phono3py,
+    mlp_params: Union[str, dict, PypolymlpParams],
+    displacement_distance: float,
+    number_of_snapshots: int,
+    random_seed: int,
+    log_level: int,
+):
+    """Run pypolymlp to compute forces."""
+    if displacement_distance is None:
+        _displacement_distance = 0.001
+    else:
+        _displacement_distance = displacement_distance
+
+    if log_level:
+        print("-" * 29 + " pypolymlp start " + "-" * 30)
+        print("MLP parameters:")
+        if mlp_params:
+            for k, v in asdict(parse_mlp_params(mlp_params)).items():
+                if v is not None:
+                    print(f"  {k}: {v}")
+
+    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(
+                "."
+            )
+        )
+    ph3py.generate_displacements(
+        distance=_displacement_distance,
+        is_plusminus=True,
+        number_of_snapshots=number_of_snapshots,
+        random_seed=random_seed,
+    )
+
+    if log_level:
+        print(
+            "  Number of supercells for computing forces: "
+            f"{ph3py.displacements.shape[0]}",
+            flush=True,
+        )
+
+    if ph3py.mlp_dataset is None:
+        msg = "mlp_dataset has to be set before calling this method."
+        raise RuntimeError(msg)
+    if ph3py.supercells_with_displacements is None:
+        raise RuntimeError("Displacements are not set. Run generate_displacements.")
+    ph3py.develop_mlp(params=mlp_params)
+    ph3py.evaluate_mlp()
+
+    if log_level:
+        print("-" * 30 + " pypolymlp end " + "-" * 31, flush=True)
+
+
 def _create_phono3py_fc2(
     phono3py: Phono3py,
     ph3py_yaml: Optional[Phono3pyYaml],

phono3py/cui/load.py

@@ -49,10 +49,9 @@ from phonopy.structure.cells import determinant
 
 from phono3py import Phono3py
 from phono3py.cui.create_force_constants import (
-    displacements_in_dataset,
     forces_in_dataset,
     parse_forces,
-    read_type2_dataset,
+    run_pypolymlp_to_compute_forces,
 )
 from phono3py.file_IO import read_fc2_from_hdf5, read_fc3_from_hdf5
 from phono3py.interface.phono3py_yaml import Phono3pyYaml
@@ -86,6 +85,8 @@ def load(
     symmetrize_fc: bool = True,
     is_mesh_symmetry: bool = True,
     is_compact_fc: bool = False,
+    use_pypolymlp: bool = False,
+    mlp_params: Optional[dict] = None,
     use_grg: bool = False,
     make_r0_average: bool = True,
     symprec: float = 1e-5,
@@ -236,6 +237,10 @@ def load(
             True: (primitive, supecell, 3, 3) False: (supercell, supecell, 3, 3)
         where 'supercell' and 'primitive' indicate number of atoms in these
         cells. Default is False.
+    use_pypolymlp : bool, optional
+        Use pypolymlp for generating force constants. Default is False.
+    mlp_params : dict, optional
+        A set of parameters used by machine learning potentials.
     use_grg : bool, optional
         Use generalized regular grid when True. Default is False.
     make_r0_average : bool, optional
@@ -337,6 +342,7 @@ def load(
         forces_fc3_filename=forces_fc3_filename,
         forces_fc2_filename=forces_fc2_filename,
         phono3py_yaml_filename=phono3py_yaml,
+        use_pypolymlp=use_pypolymlp,
         log_level=log_level,
     )
 
@@ -348,6 +354,8 @@ def load(
             fc_calculator_options=fc_calculator_options,
             symmetrize_fc=symmetrize_fc,
             is_compact_fc=is_compact_fc,
+            use_pypolymlp=use_pypolymlp,
+            mlp_params=mlp_params,
             log_level=log_level,
         )
 
@@ -370,6 +378,7 @@ def set_dataset_and_force_constants(
     forces_fc2_filename: Optional[Union[os.PathLike, Sequence]] = None,
     phono3py_yaml_filename: Optional[os.PathLike] = None,
     cutoff_pair_distance: Optional[float] = None,
+    use_pypolymlp: bool = False,
     log_level: int = 0,
 ) -> dict:
     """Set displacements, forces, and create force constants.
@@ -387,7 +396,7 @@ def set_dataset_and_force_constants(
 
     """
     read_fc = {"fc2": False, "fc3": False}
-    read_fc["fc3"] = _set_dataset_or_fc3(
+    read_fc["fc3"], dataset = _get_dataset_or_fc3(
         ph3py,
         ph3py_yaml=ph3py_yaml,
         fc3_filename=fc3_filename,
@@ -396,20 +405,20 @@ def set_dataset_and_force_constants(
         cutoff_pair_distance=cutoff_pair_distance,
         log_level=log_level,
     )
-    read_fc["fc2"] = _set_dataset_phonon_dataset_or_fc2(
+    if not read_fc["fc3"]:
+        if use_pypolymlp:
+            ph3py.mlp_dataset = dataset
+        else:
+            ph3py.dataset = dataset
+    read_fc["fc2"], phonon_dataset = _get_dataset_phonon_dataset_or_fc2(
         ph3py,
         ph3py_yaml=ph3py_yaml,
         fc2_filename=fc2_filename,
         forces_fc2_filename=forces_fc2_filename,
         log_level=log_level,
     )
-
+    if not read_fc["fc2"]:
-    # Cases that dataset is in phono3py.yaml but not forces.
+        ph3py.phonon_dataset = phonon_dataset
-    if ph3py.dataset is None:
-        if ph3py_yaml is not None and ph3py_yaml.dataset is not None:
-            ph3py.dataset = ph3py_yaml.dataset
-        if ph3py_yaml is not None and ph3py_yaml.phonon_dataset is not None:
-            ph3py.phonon_dataset = ph3py_yaml.phonon_dataset
 
     return read_fc
 
@@ -418,9 +427,14 @@ def compute_force_constants_from_datasets(
     ph3py: Phono3py,
     read_fc: dict,
     fc_calculator: Optional[str] = None,
-    fc_calculator_options: Optional[dict] = None,
+    fc_calculator_options: Optional[Union[dict, str]] = None,
     symmetrize_fc: bool = True,
     is_compact_fc: bool = True,
+    use_pypolymlp: bool = False,
+    mlp_params: Optional[Union[dict, str]] = None,
+    displacement_distance: Optional[float] = None,
+    number_of_snapshots: Optional[int] = None,
+    random_seed: Optional[int] = None,
     log_level: int = 0,
 ):
     """Compute force constants from datasets.
@@ -435,13 +449,25 @@ def compute_force_constants_from_datasets(
             fc2 : bool
 
     """
-    if not read_fc["fc3"] and ph3py.dataset:
+    if not read_fc["fc3"]:
+        if ph3py.dataset or ph3py.mlp_dataset:
+            if ph3py.mlp_dataset and use_pypolymlp:
+                run_pypolymlp_to_compute_forces(
+                    ph3py,
+                    mlp_params=mlp_params,
+                    displacement_distance=displacement_distance,
+                    number_of_snapshots=number_of_snapshots,
+                    random_seed=random_seed,
+                    log_level=log_level,
+                )
+
             ph3py.produce_fc3(
                 symmetrize_fc3r=symmetrize_fc,
                 is_compact_fc=is_compact_fc,
                 fc_calculator=fc_calculator,
                 fc_calculator_options=fc_calculator_options,
             )
+
             if log_level and symmetrize_fc and fc_calculator is None:
                 print("fc3 was symmetrized.")
 
@@ -467,7 +493,7 @@ def compute_force_constants_from_datasets(
                 print("fc2 was symmetrized.")
 
 
-def _set_dataset_or_fc3(
+def _get_dataset_or_fc3(
     ph3py: Phono3py,
     ph3py_yaml: Optional[Phono3pyYaml] = None,
     fc3_filename: Optional[os.PathLike] = None,
@@ -475,152 +501,115 @@ def _set_dataset_or_fc3(
     phono3py_yaml_filename: Optional[os.PathLike] = None,
     cutoff_pair_distance: Optional[float] = None,
     log_level: int = 0,
-) -> bool:
+) -> tuple[bool, dict]:
     p2s_map = ph3py.primitive.p2s_map
     read_fc3 = False
-    if fc3_filename is not None:
+    dataset = None
-        fc3 = read_fc3_from_hdf5(filename=fc3_filename, p2s_map=p2s_map)
+    if fc3_filename is not None or pathlib.Path("fc3.hdf5").exists():
-        _check_fc3_shape(ph3py, fc3, filename=fc3_filename)
+        if fc3_filename is None:
+            _fc3_filename = "fc3.hdf5"
+        else:
+            _fc3_filename = fc3_filename
+        fc3 = read_fc3_from_hdf5(filename=_fc3_filename, p2s_map=p2s_map)
+        _check_fc3_shape(ph3py, fc3, filename=_fc3_filename)
         ph3py.fc3 = fc3
         read_fc3 = True
         if log_level:
-            print('fc3 was read from "%s".' % fc3_filename)
+            print(f'fc3 was read from "{_fc3_filename}".')
-    elif forces_fc3_filename is not None:
-        force_filename = forces_fc3_filename
-        _set_dataset_for_fc3(
-            ph3py,
-            ph3py_yaml,
-            force_filename,
-            phono3py_yaml_filename,
-            cutoff_pair_distance,
-            log_level,
-        )
-    elif pathlib.Path("fc3.hdf5").exists():
-        fc3 = read_fc3_from_hdf5(filename="fc3.hdf5", p2s_map=p2s_map)
-        _check_fc3_shape(ph3py, fc3)
-        ph3py.fc3 = fc3
-        read_fc3 = True
-        if log_level:
-            print('fc3 was read from "fc3.hdf5".')
-    elif dataset := read_type2_dataset(
-        len(ph3py.supercell), "FORCES_FC3", log_level=log_level
-    ):  # := Assignment Expressions (Python>=3.8)
-        ph3py.dataset = dataset
-    elif ph3py_yaml is not None and ph3py_yaml.dataset is not None:
-        if pathlib.Path("FORCES_FC3").exists() and displacements_in_dataset(
-            ph3py_yaml.dataset
-        ):
-            _set_dataset_for_fc3(
-                ph3py,
-                ph3py_yaml,
-                "FORCES_FC3",
-                phono3py_yaml_filename,
-                cutoff_pair_distance,
-                log_level,
-            )
-        if forces_in_dataset(ph3py_yaml.dataset):
-            _set_dataset_for_fc3(
-                ph3py,
-                ph3py_yaml,
-                None,
-                phono3py_yaml_filename,
-                cutoff_pair_distance,
-                log_level,
-            )
-    return read_fc3
-
-
-def _set_dataset_phonon_dataset_or_fc2(
-    ph3py: Phono3py,
-    ph3py_yaml: Optional[Phono3pyYaml] = None,
-    fc2_filename: Optional[os.PathLike] = None,
-    forces_fc2_filename: Optional[Union[os.PathLike, Sequence]] = None,
-    log_level: int = 0,
-) -> bool:
-    phonon_p2s_map = ph3py.phonon_primitive.p2s_map
-    read_fc2 = False
-    if fc2_filename is not None:
-        fc2 = read_fc2_from_hdf5(filename=fc2_filename, p2s_map=phonon_p2s_map)
-        _check_fc2_shape(ph3py, fc2, filename=fc2_filename)
-        ph3py.fc2 = fc2
-        read_fc2 = True
-        if log_level:
-            print('fc2 was read from "%s".' % fc2_filename)
-    elif forces_fc2_filename is not None:
-        force_filename = forces_fc2_filename
-        _set_dataset_for_fc2(
-            ph3py,
-            ph3py_yaml,
-            force_filename,
-            "phonon_fc2",
-            log_level,
-        )
-    elif pathlib.Path("fc2.hdf5").exists():
-        fc2 = read_fc2_from_hdf5(filename="fc2.hdf5", p2s_map=phonon_p2s_map)
-        _check_fc2_shape(ph3py, fc2)
-        ph3py.fc2 = fc2
-        read_fc2 = True
-        if log_level:
-            print('fc2 was read from "fc2.hdf5".')
-    elif (
-        pathlib.Path("FORCES_FC2").exists()
-        and ph3py.phonon_supercell_matrix is not None
-    ):
-        _set_dataset_for_fc2(
-            ph3py,
-            ph3py_yaml,
-            "FORCES_FC2",
-            "phonon_fc2",
-            log_level,
-        )
-    elif (
-        ph3py_yaml is not None
-        and ph3py_yaml.phonon_dataset is not None
-        and forces_in_dataset(ph3py_yaml.phonon_dataset)
-    ):
-        _set_dataset_for_fc2(
-            ph3py,
-            ph3py_yaml,
-            None,
-            "phonon_fc2",
-            log_level,
-        )
     elif (
         ph3py_yaml is not None
         and ph3py_yaml.dataset is not None
         and forces_in_dataset(ph3py_yaml.dataset)
     ):
-        _set_dataset_for_fc2(
+        dataset = _get_dataset_for_fc3(
             ph3py,
             ph3py_yaml,
             None,
-            "fc2",
+            phono3py_yaml_filename,
+            cutoff_pair_distance,
+            log_level,
+        )
+    elif forces_fc3_filename is not None or pathlib.Path("FORCES_FC3").exists():
+        if forces_fc3_filename is None:
+            force_filename = "FORCES_FC3"
+        else:
+            force_filename = forces_fc3_filename
+        dataset = _get_dataset_for_fc3(
+            ph3py,
+            ph3py_yaml,
+            force_filename,
+            phono3py_yaml_filename,
+            cutoff_pair_distance,
+            log_level,
+        )
+        if not forces_in_dataset(dataset):
+            dataset = None
+
+    return read_fc3, dataset
+
+
+def _get_dataset_phonon_dataset_or_fc2(
+    ph3py: Phono3py,
+    ph3py_yaml: Optional[Phono3pyYaml] = None,
+    fc2_filename: Optional[os.PathLike] = None,
+    forces_fc2_filename: Optional[Union[os.PathLike, Sequence]] = None,
+    log_level: int = 0,
+) -> tuple[bool, dict, dict]:
+    phonon_p2s_map = ph3py.phonon_primitive.p2s_map
+    read_fc2 = False
+    phonon_dataset = None
+    if fc2_filename is not None or pathlib.Path("fc2.hdf5").exists():
+        if fc2_filename is None:
+            _fc2_filename = "fc2.hdf5"
+        else:
+            _fc2_filename = fc2_filename
+        fc2 = read_fc2_from_hdf5(filename=_fc2_filename, p2s_map=phonon_p2s_map)
+        _check_fc2_shape(ph3py, fc2, filename=_fc2_filename)
+        ph3py.fc2 = fc2
+        read_fc2 = True
+        if log_level:
+            print(f'fc2 was read from "{fc2_filename}".')
+    elif (
+        ph3py_yaml is not None
+        and ph3py_yaml.phonon_dataset is not None
+        and forces_in_dataset(ph3py_yaml.phonon_dataset)
+    ):
+        phonon_dataset = _get_dataset_for_fc2(
+            ph3py,
+            ph3py_yaml,
+            None,
+            "phonon_fc2",
             log_level,
         )
     elif (
-        pathlib.Path("FORCES_FC3").exists()
+        forces_fc2_filename is not None or pathlib.Path("FORCES_FC2").exists()
-        and ph3py.phonon_supercell_matrix is not None
+    ) and ph3py.phonon_supercell_matrix:
-    ):
+        if forces_fc2_filename is None:
-        # suppose fc3.hdf5 is read but fc2.hdf5 doesn't exist.
+            force_filename = forces_fc2_filename
-        _set_dataset_for_fc2(
+        else:
+            force_filename = forces_fc2_filename
+        phonon_dataset = _get_dataset_for_fc2(
             ph3py,
             ph3py_yaml,
-            "FORCES_FC3",
+            force_filename,
-            "fc2",
+            "phonon_fc2",
             log_level,
         )
-    return read_fc2
+        if not forces_in_dataset(phonon_dataset):
+            phonon_dataset = None
+
+    return read_fc2, phonon_dataset
 
 
-def _set_dataset_for_fc3(
+def _get_dataset_for_fc3(
     ph3py: Phono3py,
     ph3py_yaml: Optional[Phono3pyYaml],
     force_filename,
     phono3py_yaml_filename,
     cutoff_pair_distance,
     log_level,
-):
+) -> dict:
-    ph3py.dataset = parse_forces(
+    dataset = parse_forces(
         ph3py,
         ph3py_yaml=ph3py_yaml,
         cutoff_pair_distance=cutoff_pair_distance,
@@ -629,9 +618,10 @@ def _set_dataset_for_fc3(
         fc_type="fc3",
         log_level=log_level,
     )
+    return dataset
 
 
-def _set_dataset_for_fc2(
+def _get_dataset_for_fc2(
     ph3py: Phono3py,
     ph3py_yaml: Optional[Phono3pyYaml],
     force_filename,
@@ -645,10 +635,7 @@ def _set_dataset_for_fc2(
         fc_type=fc_type,
         log_level=log_level,
     )
-    if fc_type == "phonon_fc2":
+    return dataset
-        ph3py.phonon_dataset = dataset
-    else:
-        ph3py.dataset = dataset
 
 
 def _check_fc2_shape(ph3py: Phono3py, fc2, filename="fc2.hdf5"):

phono3py/cui/phono3py_argparse.py

@@ -471,6 +471,15 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
     parser.add_argument(
         "--mesh", nargs="+", dest="mesh_numbers", default=None, help="Mesh numbers"
     )
+    parser.add_argument(
+        "--mlp-params",
+        dest="mlp_params",
+        default=None,
+        help=(
+            "Parameters for machine learning potentials as comma separated "
+            "string with the style of key = values"
+        ),
+    )
     parser.add_argument(
         "--mv",
         "--mass-variances",
@@ -622,7 +631,7 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
         "--pypolymlp",
         dest="use_pypolymlp",
         action="store_true",
-        default=None,
+        default=False,
         help="Use pypolymlp and symfc for generating force constants",
     )
     parser.add_argument(
@@ -647,6 +656,21 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
         default=False,
         help="Print out smallest information",
     )
+    parser.add_argument(
+        "--random-seed",
+        dest="random_seed",
+        type=int,
+        default=None,
+        help="Random seed by a 32 bit unsigned integer",
+    )
+    parser.add_argument(
+        "--rd",
+        "--random-displacements",
+        dest="random_displacements",
+        type=int,
+        default=None,
+        help="Number of supercells with random displacements",
+    )
     parser.add_argument(
         "--read-collision",
         dest="read_collision",

phono3py/cui/phono3py_script.py

@@ -528,6 +528,7 @@ def store_force_constants(
             ph3py_yaml=ph3py_yaml,
             phono3py_yaml_filename=phono3py_yaml_filename,
             cutoff_pair_distance=settings.cutoff_pair_distance,
+            use_pypolymlp=settings.use_pypolymlp,
             log_level=log_level,
         )
         try:
@@ -538,6 +539,11 @@ def store_force_constants(
                 fc_calculator_options=fc_calculator_options,
                 symmetrize_fc=settings.fc_symmetry,
                 is_compact_fc=settings.is_compact_fc,
+                use_pypolymlp=settings.use_pypolymlp,
+                mlp_params=settings.mlp_params,
+                displacement_distance=settings.displacement_distance,
+                number_of_snapshots=settings.random_displacements,
+                random_seed=settings.random_seed,
                 log_level=log_level,
             )
         except ForceCalculatorRequiredError:
@@ -605,8 +611,8 @@ def _show_fc_calculator_not_found(log_level):
         print(
             "Built-in force constants calculator doesn't support the "
             "dispalcements-forces dataset. "
-            "An external force calculator, e.g., ALM (--alm), has to be used "
+            "An external force calculator, e.g., symfc (--symfc_ or ALM (--alm), "
-            "to compute force constants."
+            "has to be used to compute force constants."
         )
         print_error()
     sys.exit(1)

phono3py/cui/show_log.py

@@ -107,7 +107,7 @@ def show_phono3py_cells(phono3py: Phono3py):
         print_cell(phonon_primitive)
         print("-" * 21 + " supercell for harmonic phonon " + "-" * 22)
         print_cell(phonon_supercell, mapping=phonon_primitive.s2p_map)
-    print("-" * 76)
+    print("-" * 76, flush=True)
 
 
 def show_phono3py_force_constants_settings(settings: Phono3pySettings):

pyproject.toml

@@ -3,14 +3,14 @@ requires = ["setuptools", "wheel", "numpy"]
 
 [tool.ruff]
 line-length = 88
-select = [
+lint.select = [
     "F", # Flake8
     "B", # Black
     "I", # isort
     "E", # pycodestyle-errors
     "D", # pydocstyle
 ]
-extend-ignore = [
+lint.extend-ignore = [
     "D417",
     "D100",
 ]

test/api/test_api_phono3py.py

@@ -170,7 +170,7 @@ def test_type2_forces_energies_setter_Si(si_111_222_rd: Phono3py):
     np.testing.assert_allclose(ph3_in.phonon_forces + 1, ph3.phonon_forces)
 
 
-def test_use_pypolymlp_mgs(mgo_222rd_444rd: Phono3py):
+def test_use_pypolymlp_mgo(mgo_222rd_444rd: Phono3py):
     """Test use_pypolymlp in produce_fc3."""
     pytest.importorskip("pypolymlp")
 
@@ -193,12 +193,16 @@ def test_use_pypolymlp_mgs(mgo_222rd_444rd: Phono3py):
 
     atom_energies = {"Mg": -0.00896717, "O": -0.95743902}
     params = PypolymlpParams(gtinv_maxl=(4, 4), atom_energies=atom_energies)
-    ph3.produce_fc3(use_pypolymlp=True, mlp_params=params, fc_calculator="symfc")
+
+    ph3.generate_displacements(distance=0.001, is_plusminus=True)
+    ph3.develop_mlp(params=params)
+    ph3.evaluate_mlp()
+    ph3.produce_fc3(fc_calculator="symfc")
     ph3.produce_fc2(fc_calculator="symfc")
 
     ph3.mesh_numbers = 30
     ph3.init_phph_interaction()
     ph3.run_thermal_conductivity(temperatures=[300])
     assert (
-        pytest.approx(63.6001137, abs=1e-3) == ph3.thermal_conductivity.kappa[0, 0, 0]
+        pytest.approx(63.0018546, abs=1e-3) == ph3.thermal_conductivity.kappa[0, 0, 0]
     )

test/conftest.py

@@ -415,7 +415,6 @@ def nacl_pbe_cutoff_fc3(request) -> Phono3py:
             count += 1
     ph3.dataset = dataset
     ph3.produce_fc3()
-    # ph3.produce_fc3(symmetrize_fc3r=True)
     return ph3
 
 
@@ -504,7 +503,7 @@ def mgo_222rd_444rd() -> Phono3py:
     4 and 400 supercells for fc2 and fc3, respectively.
 
     """
-    yaml_filename = cwd / "phono3py_params_MgO-222rd-444fd.yaml.xz"
+    yaml_filename = cwd / "phono3py_params_MgO-222rd-444rd.yaml.xz"
     return phono3py.load(yaml_filename, produce_fc=False, log_level=1)