Write fc3_nonzero_indices in fc3.hdf5

phono3py/api_phono3py.py

@@ -338,7 +338,7 @@ class Phono3py:
     def fc3(self) -> NDArray | None:
         """Setter and getter of third order force constants (fc3).
 
-        ndarray
+        ndarray, optional
             fc3 shape is either (supercell, supercell, supercell, 3, 3, 3) or
             (primitive, supercell, supercell, 3, 3, 3),
             where 'supercell' and 'primitive' indicate number of atoms in

phono3py/cui/create_force_constants.py

@@ -142,6 +142,7 @@ def create_phono3py_force_constants(
                     fc_calculator_options, 3
                 ),
             )
+        assert phono3py.fc3 is not None, "fc3 is not set."
 
         cutoff_distance = settings.cutoff_fc3_distance
         if cutoff_distance is not None and cutoff_distance > 0:
@@ -157,9 +158,10 @@ def create_phono3py_force_constants(
             else:
                 filename = "fc3." + output_filename + ".hdf5"
             if log_level:
-                print('Writing fc3 to "%s".' % filename)
+                print(f'Writing fc3 to "{filename}".')
             write_fc3_to_hdf5(
                 phono3py.fc3,
+                fc3_nonzero_indices=phono3py.fc3_nonzero_indices,
                 filename=filename,
                 p2s_map=phono3py.primitive.p2s_map,
                 compression=settings.hdf5_compression,

phono3py/cui/phono3py_script.py

@@ -679,6 +679,7 @@ def _store_force_constants(ph3py: Phono3py, settings: Phono3pySettings, log_leve
     if not read_fc3:
         write_fc3_to_hdf5(
             ph3py.fc3,
+            fc3_nonzero_indices=ph3py.fc3_nonzero_indices,
             p2s_map=ph3py.primitive.p2s_map,
             compression=settings.hdf5_compression,
         )

phono3py/file_IO.py

@@ -42,6 +42,7 @@ from typing import Optional, TextIO, Union
 
 import h5py
 import numpy as np
+from numpy.typing import NDArray
 from phonopy.cui.load_helper import read_force_constants_from_hdf5
 
 # This import is deactivated for a while.
@@ -300,29 +301,42 @@ def _write_FORCES_FC3_typeI(
             count += 1
 
 
-def write_fc3_to_hdf5(fc3, filename="fc3.hdf5", p2s_map=None, compression="gzip"):
+def write_fc3_to_hdf5(
+    fc3: NDArray,
+    fc3_nonzero_indices: NDArray | None = None,
+    filename: str = "fc3.hdf5",
+    p2s_map: NDArray | None = None,
+    compression: str = "gzip",
+):
     """Write fc3 in fc3.hdf5.
 
     Parameters
     ----------
-    force_constants : ndarray
-        Force constants
-        shape=(n_satom, n_satom, n_satom, 3, 3, 3) or
-        (n_patom, n_satom, n_satom,3,3,3), dtype=double
+    fc3 : ndarray
+        Force constants shape=(n_satom, n_satom, n_satom, 3, 3, 3) or (n_patom,
+        n_satom, n_satom,3,3,3), dtype=double
+    fc3_nonzero_indices : ndarray, optional
+        Non-zero indices of fc3. shape=(n_satom, n_satom, n_satom) or (n_patom,
+        n_satom, n_satom), dtype="byte". If this is given, fc3 is in compact
+        format. Otherwise, it is in full format.
     filename : str
         Filename to be used.
     p2s_map : ndarray, optional
-        Primitive atom indices in supercell index system
-        shape=(n_patom,), dtype=intc
+        Primitive atom indices in supercell index system shape=(n_patom,),
+        dtype=intc
     compression : str or int, optional
-        h5py's lossless compression filters (e.g., "gzip", "lzf"). None gives
-        no compression. See the detail at docstring of
-        h5py.Group.create_dataset. Default is "gzip".
+        h5py's lossless compression filters (e.g., "gzip", "lzf"). None gives no
+        compression. See the detail at docstring of h5py.Group.create_dataset.
+        Default is "gzip".
 
     """
     with h5py.File(filename, "w") as w:
         w.create_dataset("version", data=np.bytes_(__version__))
         w.create_dataset("fc3", data=fc3, compression=compression)
+        if fc3_nonzero_indices is not None:
+            w.create_dataset(
+                "fc3_nonzero_indices", data=fc3_nonzero_indices, compression=compression
+            )
         if p2s_map is not None:
             w.create_dataset("p2s_map", data=p2s_map)
 

test/cui/test_phono3py_load_script.py

@@ -71,8 +71,13 @@ def test_phono3py_load():
             file_path.unlink()
 
 
-@pytest.mark.parametrize("fc_calculator,exit_code", [(None, 0), ("symfc", 0)])
-def test_phono3py_load_with_typeII_dataset(fc_calculator, exit_code):
+@pytest.mark.parametrize(
+    "fc_calculator,fc_calculator_options",
+    [(None, None), ("symfc", None), ("symfc", "|cutoff=4.0")],
+)
+def test_phono3py_load_with_typeII_dataset(
+    fc_calculator: str | None, fc_calculator_options: str | None
+):
     """Test phono3py-load script with typeII dataset.
 
     When None, fallback to symfc.
@@ -80,16 +85,24 @@ def test_phono3py_load_with_typeII_dataset(fc_calculator, exit_code):
     """
     pytest.importorskip("symfc")
     argparse_control = _get_phono3py_load_args(
-        cwd / ".." / "phono3py_params-Si111-rd.yaml.xz", fc_calculator=fc_calculator
+        cwd / ".." / "phono3py_params-Si111-rd.yaml.xz",
+        fc_calculator=fc_calculator,
+        fc_calculator_options=fc_calculator_options,
     )
     with pytest.raises(SystemExit) as excinfo:
         main(**argparse_control)
-    assert excinfo.value.code == exit_code
+    assert excinfo.value.code == 0
 
     # Clean files created by phono3py-load script.
     for created_filename in ("phono3py.yaml", "fc2.hdf5", "fc3.hdf5"):
         file_path = pathlib.Path(cwd_called / created_filename)
         if file_path.exists():
+            if created_filename == "fc3.hdf5":
+                with h5py.File(file_path, "r") as f:
+                    if fc_calculator_options is None:
+                        assert "fc3_nonzero_indices" not in f
+                    else:
+                        assert "fc3_nonzero_indices" in f
             file_path.unlink()