Merge pull request #254 from phonopy/fix-gen-disps-wo-forces

Fix phono3py.load: loading dataset wo forces

.github/workflows/phono3py-pytest-conda.yml

@@ -12,7 +12,7 @@ jobs:
         shell: bash -l {0}
     strategy:
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
 
     steps:
     - uses: actions/checkout@v4

.github/workflows/publish-to-test-pypi.yml

@@ -11,7 +11,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.12, ]
+        python-version: ["3.12"]
 
     steps:
     - uses: actions/checkout@v4

phono3py/api_phono3py.py

@@ -81,7 +81,7 @@
 from phono3py.interface.fc_calculator import get_fc3
 from phono3py.interface.phono3py_yaml import Phono3pyYaml
 from phono3py.phonon.grid import BZGrid
-from phono3py.phonon3.dataset import get_displacements_and_forces_fc3
+from phono3py.phonon3.dataset import forces_in_dataset, get_displacements_and_forces_fc3
 from phono3py.phonon3.displacement_fc3 import (
     direction_to_displacement,
     get_third_order_displacements,
@@ -846,6 +846,9 @@ def displacements(self):
         """
         dataset = self._dataset
 
+        if self._dataset is None:
+            raise RuntimeError("displacement dataset is not set.")
+
         if "first_atoms" in dataset:
             num_scells = len(dataset["first_atoms"])
             for disp1 in dataset["first_atoms"]:
@@ -863,7 +866,7 @@ def displacements(self):
                 for disp2 in disp1["second_atoms"]:
                     displacements[i, disp2["number"]] = disp2["displacement"]
                     i += 1
-        elif "forces" in dataset or "displacements" in dataset:
+        elif "displacements" in dataset:
             displacements = dataset["displacements"]
         else:
             raise RuntimeError("displacement dataset has wrong format.")
@@ -1218,10 +1221,10 @@ def run_phonon_solver(self, grid_points=None):
 
     def generate_displacements(
         self,
-        distance=0.03,
-        cutoff_pair_distance=None,
-        is_plusminus="auto",
-        is_diagonal=True,
+        distance: float = 0.03,
+        cutoff_pair_distance: Optional[float] = None,
+        is_plusminus: Union[bool, str] = "auto",
+        is_diagonal: bool = True,
         number_of_snapshots: Optional[int] = None,
         random_seed: Optional[int] = None,
         is_random_distance: bool = False,
@@ -1410,7 +1413,7 @@ def produce_fc3(
         self,
         symmetrize_fc3r: bool = False,
         is_compact_fc: bool = False,
-        fc_calculator: Optional[str] = None,
+        fc_calculator: Optional[Union[str, dict]] = None,
         fc_calculator_options: Optional[Union[str, dict]] = None,
     ):
         """Calculate fc3 from displacements and forces.
@@ -2150,7 +2153,9 @@ def run_thermal_conductivity(
                 log_level=_log_level,
             )
 
-    def save(self, filename="phono3py_params.yaml", settings=None):
+    def save(
+        self, filename: str = "phono3py_params.yaml", settings: Optional[dict] = None
+    ):
         """Save parameters in Phono3py instants into file.
 
         Parameters
@@ -2541,6 +2546,11 @@ def _get_forces_energies(
         Return None if tagert data is not found rather than raising exception.
 
         """
+        if self._dataset is None:
+            return None
+        if not forces_in_dataset(self._dataset):
+            return None
+
         if target in self._dataset:  # type-2
             return self._dataset[target]
         elif "first_atoms" in self._dataset:  # type-1

phono3py/cui/create_force_constants.py

@@ -68,6 +68,7 @@
 )
 from phono3py.interface.fc_calculator import extract_fc2_fc3_calculators
 from phono3py.interface.phono3py_yaml import Phono3pyYaml
+from phono3py.phonon3.dataset import forces_in_dataset
 from phono3py.phonon3.fc3 import (
     set_permutation_symmetry_fc3,
     set_translational_invariance_fc3,
@@ -95,8 +96,6 @@ def create_phono3py_force_constants(
     symmetrize_fc3r = settings.is_symmetrize_fc3_r or settings.fc_symmetry
     symmetrize_fc2 = settings.is_symmetrize_fc2 or settings.fc_symmetry
 
-    (fc_calculator, fc_calculator_options) = get_fc_calculator_params(settings)
-
     if log_level:
         show_phono3py_force_constants_settings(settings)
 
@@ -113,6 +112,9 @@ def create_phono3py_force_constants(
     ):
         pass
     else:
+        (fc_calculator, fc_calculator_options) = get_fc_calculator_params(
+            settings, log_level=(not settings.read_fc3) * 1
+        )
         if settings.read_fc3:
             _read_phono3py_fc3(phono3py, symmetrize_fc3r, input_filename, log_level)
         else:  # fc3 from FORCES_FC3 or ph3py_yaml
@@ -292,10 +294,17 @@ def parse_forces(
     # Type-1 FORCES_FC*.
     # dataset comes either from disp_fc*.yaml or phono3py*.yaml.
     if not forces_in_dataset(dataset):
-        if fc_type == "phonon_fc2":
-            parse_FORCES_FC2(dataset, filename=force_filename)
-        else:
-            parse_FORCES_FC3(dataset, filename=force_filename)
+        if force_filename is not None:
+            if fc_type == "phonon_fc2":
+                parse_FORCES_FC2(dataset, filename=force_filename)
+            else:
+                parse_FORCES_FC3(dataset, filename=force_filename)
+
+            if log_level:
+                print(
+                    f'Sets of supercell forces were read from "{force_filename}".',
+                    flush=True,
+                )
 
         # Unit of displacements is already converted.
         # Therefore, only unit of forces is converted.
@@ -305,28 +314,10 @@ def parse_forces(
                 force_to_eVperA=physical_units["force_to_eVperA"],
             )
 
-        if log_level:
-            print('Sets of supercell forces were read from "%s".' % force_filename)
-            sys.stdout.flush()
-
     return dataset
 
 
-def forces_in_dataset(dataset: dict) -> bool:
-    """Return whether forces in dataset or not."""
-    return "forces" in dataset or (
-        "first_atoms" in dataset and "forces" in dataset["first_atoms"][0]
-    )
-
-
-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
-
-
-def get_fc_calculator_params(settings):
+def get_fc_calculator_params(settings, log_level=0):
     """Return fc_calculator and fc_calculator_params from settings."""
     fc_calculator = None
     fc_calculator_list = []
@@ -339,13 +330,49 @@ def get_fc_calculator_params(settings):
         if fc_calculator_list:
             fc_calculator = "|".join(fc_calculator_list)
 
-    fc_calculator_options = None
-    if settings.fc_calculator_options is not None:
-        fc_calculator_options = settings.fc_calculator_options
+    fc_calculator_options = settings.fc_calculator_options
+    if settings.cutoff_pair_distance:
+        if fc_calculator_list and fc_calculator_list[-1] in ("alm", "symfc"):
+            if fc_calculator_list[-1] == "alm":
+                cutoff_str = f"-1 {settings.cutoff_pair_distance}"
+            if fc_calculator_list[-1] == "symfc":
+                cutoff_str = f"{settings.cutoff_pair_distance}"
+            fc_calculator_options = _set_cutoff_in_fc_calculator_options(
+                fc_calculator_options,
+                cutoff_str,
+                log_level,
+            )
 
     return fc_calculator, fc_calculator_options
 
 
+def _set_cutoff_in_fc_calculator_options(
+    fc_calculator_options: Optional[str],
+    cutoff_str: str,
+    log_level: int,
+):
+    str_appended = f"cutoff={cutoff_str}"
+    calc_opts = fc_calculator_options
+    if calc_opts is None:
+        calc_opts = "|"
+    if "|" in calc_opts:
+        calc_opts_fc2, calc_opts_fc3 = [v.strip() for v in calc_opts.split("|")][:2]
+    else:
+        calc_opts_fc2 = calc_opts
+        calc_opts_fc3 = calc_opts
+
+    if calc_opts_fc3 == "":
+        calc_opts_fc3 += f"{str_appended}"
+        if log_level:
+            print(f'Set "{str_appended}" to fc_calculator_options for fc3.')
+    elif "cutoff" not in calc_opts_fc3:
+        calc_opts_fc3 += f", {str_appended}"
+        if log_level:
+            print(f'Appended "{str_appended}" to fc_calculator_options for fc3.')
+
+    return f"{calc_opts_fc2}|{calc_opts_fc3}"
+
+
 def _read_phono3py_fc3(phono3py: Phono3py, symmetrize_fc3r, input_filename, log_level):
     if input_filename is None:
         filename = "fc3.hdf5"
@@ -417,7 +444,7 @@ def _read_phono3py_fc2(phono3py, symmetrize_fc2, input_filename, log_level):
 
 def read_type2_dataset(natom, filename="FORCES_FC3", log_level=0) -> Optional[dict]:
     """Read type-2 FORCES_FC3."""
-    if not pathlib.Path(filename).exists():
+    if filename is None or not pathlib.Path(filename).exists():
         return None
 
     with open(filename, "r") as f:

phono3py/cui/load.py

@@ -49,13 +49,13 @@
 
 from phono3py import Phono3py
 from phono3py.cui.create_force_constants import (
-    forces_in_dataset,
     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
 from phono3py.interface.phono3py_yaml import Phono3pyYaml
+from phono3py.phonon3.dataset import forces_in_dataset
 from phono3py.phonon3.fc3 import show_drift_fc3
 
 
@@ -453,7 +453,7 @@ def compute_force_constants_from_datasets(
     fc3_calculator = extract_fc2_fc3_calculators(fc_calculator, 3)
     fc2_calculator = extract_fc2_fc3_calculators(fc_calculator, 2)
     if not read_fc["fc3"] and (ph3py.dataset or ph3py.mlp_dataset):
-        if ph3py.mlp_dataset and use_pypolymlp:
+        if use_pypolymlp and forces_in_dataset(ph3py.mlp_dataset):
             run_pypolymlp_to_compute_forces(
                 ph3py,
                 mlp_params=mlp_params,
@@ -462,17 +462,22 @@ def compute_force_constants_from_datasets(
                 random_seed=random_seed,
                 log_level=log_level,
             )
-        ph3py.produce_fc3(
-            symmetrize_fc3r=symmetrize_fc,
-            is_compact_fc=is_compact_fc,
-            fc_calculator=fc3_calculator,
-            fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 3),
-        )
+        if forces_in_dataset(ph3py.dataset):
+            ph3py.produce_fc3(
+                symmetrize_fc3r=symmetrize_fc,
+                is_compact_fc=is_compact_fc,
+                fc_calculator=fc3_calculator,
+                fc_calculator_options=extract_fc2_fc3_calculators(
+                    fc_calculator_options, 3
+                ),
+            )
 
-        if log_level and symmetrize_fc and fc_calculator is None:
-            print("fc3 was symmetrized.")
+            if log_level and symmetrize_fc and fc_calculator is None:
+                print("fc3 was symmetrized.")
 
-    if not read_fc["fc2"] and (ph3py.dataset or ph3py.phonon_dataset):
+    if not read_fc["fc2"] and (
+        forces_in_dataset(ph3py.dataset) or forces_in_dataset(ph3py.phonon_dataset)
+    ):
         ph3py.produce_fc2(
             symmetrize_fc2=symmetrize_fc,
             is_compact_fc=is_compact_fc,
@@ -495,6 +500,7 @@ def _get_dataset_or_fc3(
     p2s_map = ph3py.primitive.p2s_map
     read_fc3 = False
     dataset = None
+
     if fc3_filename is not None or pathlib.Path("fc3.hdf5").exists():
         if fc3_filename is None:
             _fc3_filename = "fc3.hdf5"
@@ -532,8 +538,17 @@ def _get_dataset_or_fc3(
             cutoff_pair_distance,
             log_level,
         )
-        if not forces_in_dataset(dataset):
-            dataset = None
+    elif ph3py_yaml is not None and ph3py_yaml.dataset is not None:
+        # not forces_in_dataset(ph3py_yaml.dataset)
+        # but want to read displacement dataset.
+        dataset = _get_dataset_for_fc3(
+            ph3py,
+            ph3py_yaml,
+            None,
+            phono3py_yaml_filename,
+            cutoff_pair_distance,
+            log_level,
+        )
 
     return read_fc3, dataset
 
@@ -585,8 +600,16 @@ def _get_dataset_phonon_dataset_or_fc2(
             "phonon_fc2",
             log_level,
         )
-        if not forces_in_dataset(phonon_dataset):
-            phonon_dataset = None
+    elif ph3py_yaml is not None and ph3py_yaml.phonon_dataset is not None:
+        # not forces_in_dataset(ph3py_yaml.dataset)
+        # but want to read displacement dataset.
+        phonon_dataset = _get_dataset_for_fc2(
+            ph3py,
+            ph3py_yaml,
+            None,
+            "phonon_fc2",
+            log_level,
+        )
 
     return read_fc2, phonon_dataset
 

phono3py/cui/phono3py_script.py

@@ -528,8 +528,6 @@ def store_force_constants(
         if log_level:
             print("-" * 29 + " Force constants " + "-" * 30)
 
-        (fc_calculator, fc_calculator_options) = get_fc_calculator_params(settings)
-
         read_fc = set_dataset_and_force_constants(
             phono3py,
             ph3py_yaml=ph3py_yaml,
@@ -538,6 +536,9 @@ def store_force_constants(
             use_pypolymlp=settings.use_pypolymlp,
             log_level=log_level,
         )
+        (fc_calculator, fc_calculator_options) = get_fc_calculator_params(
+            settings, log_level=(not read_fc["fc3"]) * 1
+        )
         try:
             compute_force_constants_from_datasets(
                 phono3py,

phono3py/phonon3/dataset.py

@@ -34,6 +34,8 @@
 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
 
+from typing import Optional
+
 import numpy as np
 
 
@@ -94,3 +96,12 @@ def get_displacements_and_forces_fc3(disp_dataset):
         return disp_dataset["displacements"], disp_dataset["forces"]
     else:
         raise RuntimeError("disp_dataset doesn't contain correct information.")
+
+
+def forces_in_dataset(dataset: Optional[dict]) -> bool:
+    """Return whether forces in dataset or not."""
+    if dataset is None:
+        return False
+    return "forces" in dataset or (
+        "first_atoms" in dataset and "forces" in dataset["first_atoms"][0]
+    )