Enabled cutoff-fc3-zero with comapct fc

phono3py/api_phono3py.py

@@ -1717,6 +1717,12 @@ class Phono3py:
     def cutoff_fc3_by_zero(self, cutoff_distance, fc3=None):
         """Set zero to fc3 elements out of cutoff distance.
 
+        Note
+        ----
+        fc3 is overwritten.
+
+        Parameters
+        ----------
         cutoff_distance : float
             After creating force constants, fc elements where any pair
             distance in atom triplets larger than cutoff_distance are set zero.
@@ -1727,7 +1733,11 @@ class Phono3py:
         else:
             _fc3 = fc3
         cutoff_fc3_by_zero(
-            _fc3, self._supercell, cutoff_distance, self._symprec  # overwritten
+            _fc3,
+            self._supercell,
+            cutoff_distance,
+            p2s_map=self._primitive.p2s_map,
+            symprec=self._symprec,
         )
 
     def set_permutation_symmetry(self):

phono3py/cui/create_force_constants.py

@@ -46,6 +46,7 @@ from phonopy.harmonic.force_constants import (
 )
 from phonopy.interface.calculator import get_default_physical_units
 
+from phono3py import Phono3py
 from phono3py.cui.show_log import show_phono3py_force_constants_settings
 from phono3py.file_IO import (
     get_length_of_first_line,
@@ -67,7 +68,7 @@ from phono3py.phonon3.fc3 import (
 
 
 def create_phono3py_force_constants(
-    phono3py,
+    phono3py: Phono3py,
     settings,
     ph3py_yaml=None,
     input_filename=None,
@@ -112,6 +113,16 @@ def create_phono3py_force_constants(
                 settings.fc_calculator_options,
                 log_level,
             )
+
+        cutoff_distance = settings.cutoff_fc3_distance
+        if cutoff_distance is not None and cutoff_distance > 0:
+            if log_level:
+                print(
+                    "Cutting-off fc3 by zero (cut-off distance: %f)" % cutoff_distance
+                )
+            phono3py.cutoff_fc3_by_zero(cutoff_distance)
+
+        if not settings.read_fc3:
             if output_filename is None:
                 filename = "fc3.hdf5"
             else:
@@ -125,14 +136,6 @@ def create_phono3py_force_constants(
                 compression=settings.hdf5_compression,
             )
 
-        cutoff_distance = settings.cutoff_fc3_distance
-        if cutoff_distance is not None and cutoff_distance > 0:
-            if log_level:
-                print(
-                    "Cutting-off fc3 by zero (cut-off distance: %f)" % cutoff_distance
-                )
-            phono3py.cutoff_fc3_by_zero(cutoff_distance)
-
         if log_level:
             show_drift_fc3(phono3py.fc3, primitive=phono3py.primitive)
 

phono3py/cui/phono3py_script.py

@@ -742,7 +742,7 @@ def grid_addresses_to_grid_points(grid_addresses, bz_grid):
 
 
 def store_force_constants(
-    phono3py,
+    phono3py: Phono3py,
     settings,
     ph3py_yaml,
     input_filename,
@@ -777,6 +777,14 @@ def store_force_constants(
             print_error()
             sys.exit(1)
 
+        cutoff_distance = settings.cutoff_fc3_distance
+        if cutoff_distance is not None and cutoff_distance > 0:
+            if log_level:
+                print(
+                    "Cutting-off fc3 by zero (cut-off distance: %f)" % cutoff_distance
+                )
+            phono3py.cutoff_fc3_by_zero(cutoff_distance)
+
         if not read_fc["fc3"]:
             write_fc3_to_hdf5(
                 phono3py.fc3,

phono3py/phonon3/fc3.py

@@ -499,7 +499,7 @@ def _cutoff_fc3_for_cutoff_pairs(fc3, supercell, disp_dataset, symmetry, verbose
                 _copy_permutation_symmetry_fc3_elem(fc3, ave_fc3, i, j, k)
 
 
-def cutoff_fc3_by_zero(fc3, supercell, cutoff_distance, symprec=1e-5):
+def cutoff_fc3_by_zero(fc3, supercell, cutoff_distance, p2s_map=None, symprec=1e-5):
     """Set zero in fc3 elements where pair distances are larger than cutoff."""
     num_atom = len(supercell)
     lattice = supercell.cell.T
@@ -512,11 +512,19 @@ def cutoff_fc3_by_zero(fc3, supercell, cutoff_distance, symprec=1e-5):
                 )
             )
 
-    for i, j, k in np.ndindex(num_atom, num_atom, num_atom):
-        for pair in ((i, j), (j, k), (k, i)):
-            if min_distances[pair] > cutoff_distance:
-                fc3[i, j, k] = 0
-                break
+    if fc3.shape[0] == fc3.shape[1]:
+        _p2s_map = np.arange(num_atom)
+    elif p2s_map is None or len(p2s_map) != fc3.shape[0]:
+        raise RuntimeError("Array shape of fc3 is incorrect.")
+    else:
+        _p2s_map = p2s_map
+
+    for i_index, i in enumerate(_p2s_map):
+        for j, k in np.ndindex(num_atom, num_atom):
+            for pair in ((i, j), (j, k), (k, i)):
+                if min_distances[pair] > cutoff_distance:
+                    fc3[i_index, j, k] = 0
+                    break
 
 
 def show_drift_fc3(fc3, primitive=None, name="fc3"):

test/conftest.py

@@ -192,6 +192,24 @@ def si_pbesol_iterha_111():
 def nacl_pbe(request):
     """Return Phono3py instance of NaCl 2x2x2.
 
+    * with symmetry
+    * full fc
+
+    """
+    yaml_filename = os.path.join(current_dir, "phono3py_params_NaCl222.yaml.xz")
+    enable_v2 = request.config.getoption("--v1")
+    return phono3py.load(
+        yaml_filename,
+        store_dense_gp_map=enable_v2,
+        store_dense_svecs=enable_v2,
+        log_level=1,
+    )
+
+
+@pytest.fixture(scope="session")
+def nacl_pbe_compact_fc(request):
+    """Return Phono3py instance of NaCl 2x2x2.
+
     * with symmetry
     * compact fc
 
@@ -202,6 +220,7 @@ def nacl_pbe(request):
         yaml_filename,
         store_dense_gp_map=enable_v2,
         store_dense_svecs=enable_v2,
+        is_compact_fc=True,
         log_level=1,
     )
 

test/phonon3/test_fc3.py

@@ -21,6 +21,15 @@ def test_cutoff_fc3_zero(nacl_pbe):
     assert np.abs(5259.2234163 - abs_delta) < 1e-3
 
 
+def test_cutoff_fc3_zero_compact_fc(nacl_pbe_compact_fc):
+    """Test for abrupt cut of fc3 by distance."""
+    ph = nacl_pbe_compact_fc
+    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
+
+
 def test_fc3(si_pbesol_111):
     """Test fc3 with Si PBEsol 1x1x1."""
     ph = si_pbesol_111