diff --git a/c/_phono3py.cpp b/c/_phono3py.cpp
index 139c9a76..032d2393 100644
--- a/c/_phono3py.cpp
+++ b/c/_phono3py.cpp
@@ -50,18 +50,17 @@ static Darray *convert_to_darray(nb::ndarray<> npyary) {
 //     printf("Data shift:%lu [%lu, %lu]\n", adrs_shift, i_sigma, i_temp);
 // }
 
-void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
-                        nb::ndarray<> py_g_zero, nb::ndarray<> py_frequencies,
-                        nb::ndarray<> py_eigenvectors,
-                        nb::ndarray<> py_triplets,
-                        nb::ndarray<> py_bz_grid_addresses,
-                        nb::ndarray<> py_D_diag, nb::ndarray<> py_Q,
-                        nb::ndarray<> py_fc3, nb::ndarray<> py_svecs,
-                        nb::ndarray<> py_multi, nb::ndarray<> py_masses,
-                        nb::ndarray<> py_p2s_map, nb::ndarray<> py_s2p_map,
-                        nb::ndarray<> py_band_indices, int64_t symmetrize_fc3_q,
-                        int64_t make_r0_average, nb::ndarray<> py_all_shortest,
-                        double cutoff_frequency, int64_t openmp_per_triplets) {
+void py_get_interaction(
+    nb::ndarray<> py_fc3_normal_squared, nb::ndarray<> py_g_zero,
+    nb::ndarray<> py_frequencies, nb::ndarray<> py_eigenvectors,
+    nb::ndarray<> py_triplets, nb::ndarray<> py_bz_grid_addresses,
+    nb::ndarray<> py_D_diag, nb::ndarray<> py_Q, nb::ndarray<> py_fc3,
+    nb::ndarray<> py_fc3_nonzero_indices, nb::ndarray<> py_svecs,
+    nb::ndarray<> py_multi, nb::ndarray<> py_masses, nb::ndarray<> py_p2s_map,
+    nb::ndarray<> py_s2p_map, nb::ndarray<> py_band_indices,
+    int64_t symmetrize_fc3_q, int64_t make_r0_average,
+    nb::ndarray<> py_all_shortest, double cutoff_frequency,
+    int64_t openmp_per_triplets) {
     Darray *fc3_normal_squared;
     Darray *freqs;
     _lapack_complex_double *eigvecs;
@@ -72,6 +71,7 @@ void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
     int64_t *D_diag;
     int64_t (*Q)[3];
     double *fc3;
+    char *fc3_nonzero_indices;
     double (*svecs)[3];
     int64_t (*multi)[2];
     double *masses;
@@ -100,6 +100,7 @@ void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
     } else {
         is_compact_fc3 = 1;
     }
+    fc3_nonzero_indices = (char *)py_fc3_nonzero_indices.data();
     svecs = (double (*)[3])py_svecs.data();
     for (i = 0; i < 2; i++) {
         multi_dims[i] = py_multi.shape(i);
@@ -113,9 +114,10 @@ void py_get_interaction(nb::ndarray<> py_fc3_normal_squared,
 
     ph3py_get_interaction(fc3_normal_squared, g_zero, freqs, eigvecs, triplets,
                           num_triplets, bz_grid_addresses, D_diag, Q, fc3,
-                          is_compact_fc3, svecs, multi_dims, multi, masses, p2s,
-                          s2p, band_indices, symmetrize_fc3_q, make_r0_average,
-                          all_shortest, cutoff_frequency, openmp_per_triplets);
+                          fc3_nonzero_indices, is_compact_fc3, svecs,
+                          multi_dims, multi, masses, p2s, s2p, band_indices,
+                          symmetrize_fc3_q, make_r0_average, all_shortest,
+                          cutoff_frequency, openmp_per_triplets);
 
     free(fc3_normal_squared);
     fc3_normal_squared = NULL;
@@ -129,8 +131,9 @@ void py_get_pp_collision(
     nb::ndarray<> py_triplets, nb::ndarray<> py_triplet_weights,
     nb::ndarray<> py_bz_grid_addresses, nb::ndarray<> py_bz_map,
     int64_t bz_grid_type, nb::ndarray<> py_D_diag, nb::ndarray<> py_Q,
-    nb::ndarray<> py_fc3, nb::ndarray<> py_svecs, nb::ndarray<> py_multi,
-    nb::ndarray<> py_masses, nb::ndarray<> py_p2s_map, nb::ndarray<> py_s2p_map,
+    nb::ndarray<> py_fc3, nb::ndarray<> py_fc3_nonzero_indices,
+    nb::ndarray<> py_svecs, nb::ndarray<> py_multi, nb::ndarray<> py_masses,
+    nb::ndarray<> py_p2s_map, nb::ndarray<> py_s2p_map,
     nb::ndarray<> py_band_indices, nb::ndarray<> py_temperatures_THz,
     int64_t is_NU, int64_t symmetrize_fc3_q, int64_t make_r0_average,
     nb::ndarray<> py_all_shortest, double cutoff_frequency,
@@ -147,6 +150,7 @@ void py_get_pp_collision(
     int64_t *D_diag;
     int64_t (*Q)[3];
     double *fc3;
+    char *fc3_nonzero_indices;
     double (*svecs)[3];
     int64_t (*multi)[2];
     double *masses;
@@ -176,6 +180,7 @@ void py_get_pp_collision(
     } else {
         is_compact_fc3 = 1;
     }
+    fc3_nonzero_indices = (char *)py_fc3_nonzero_indices.data();
     svecs = (double (*)[3])py_svecs.data();
     for (i = 0; i < 2; i++) {
         multi_dims[i] = py_multi.shape(i);
@@ -191,9 +196,10 @@ void py_get_pp_collision(
     ph3py_get_pp_collision(
         gamma, relative_grid_address, frequencies, eigenvectors, triplets,
         num_triplets, triplet_weights, bz_grid_addresses, bz_map, bz_grid_type,
-        D_diag, Q, fc3, is_compact_fc3, svecs, multi_dims, multi, masses, p2s,
-        s2p, band_indices, temperatures_THz, is_NU, symmetrize_fc3_q,
-        make_r0_average, all_shortest, cutoff_frequency, openmp_per_triplets);
+        D_diag, Q, fc3, fc3_nonzero_indices, is_compact_fc3, svecs, multi_dims,
+        multi, masses, p2s, s2p, band_indices, temperatures_THz, is_NU,
+        symmetrize_fc3_q, make_r0_average, all_shortest, cutoff_frequency,
+        openmp_per_triplets);
 
     free(band_indices);
     band_indices = NULL;
@@ -206,7 +212,8 @@ void py_get_pp_collision_with_sigma(
     nb::ndarray<> py_frequencies, nb::ndarray<> py_eigenvectors,
     nb::ndarray<> py_triplets, nb::ndarray<> py_triplet_weights,
     nb::ndarray<> py_bz_grid_addresses, nb::ndarray<> py_D_diag,
-    nb::ndarray<> py_Q, nb::ndarray<> py_fc3, nb::ndarray<> py_svecs,
+    nb::ndarray<> py_Q, nb::ndarray<> py_fc3,
+    nb::ndarray<> py_fc3_nonzero_indices, nb::ndarray<> py_svecs,
     nb::ndarray<> py_multi, nb::ndarray<> py_masses, nb::ndarray<> py_p2s_map,
     nb::ndarray<> py_s2p_map, nb::ndarray<> py_band_indices,
     nb::ndarray<> py_temperatures_THz, int64_t is_NU, int64_t symmetrize_fc3_q,
@@ -222,6 +229,7 @@ void py_get_pp_collision_with_sigma(
     int64_t *D_diag;
     int64_t (*Q)[3];
     double *fc3;
+    char *fc3_nonzero_indices;
     double (*svecs)[3];
     int64_t (*multi)[2];
     double *masses;
@@ -249,6 +257,7 @@ void py_get_pp_collision_with_sigma(
     } else {
         is_compact_fc3 = 1;
     }
+    fc3_nonzero_indices = (char *)py_fc3_nonzero_indices.data();
     svecs = (double (*)[3])py_svecs.data();
     for (i = 0; i < 2; i++) {
         multi_dims[i] = py_multi.shape(i);
@@ -264,8 +273,8 @@ void py_get_pp_collision_with_sigma(
     ph3py_get_pp_collision_with_sigma(
         gamma, sigma, sigma_cutoff, frequencies, eigenvectors, triplets,
         num_triplets, triplet_weights, bz_grid_addresses, D_diag, Q, fc3,
-        is_compact_fc3, svecs, multi_dims, multi, masses, p2s, s2p,
-        band_indices, temperatures_THz, is_NU, symmetrize_fc3_q,
+        fc3_nonzero_indices, is_compact_fc3, svecs, multi_dims, multi, masses,
+        p2s, s2p, band_indices, temperatures_THz, is_NU, symmetrize_fc3_q,
         make_r0_average, all_shortest, cutoff_frequency, openmp_per_triplets);
 
     free(band_indices);
diff --git a/c/phono3py.c b/c/phono3py.c
index b3793a99..47029370 100644
--- a/c/phono3py.c
+++ b/c/phono3py.c
@@ -63,12 +63,13 @@ int64_t ph3py_get_interaction(
     const _lapack_complex_double *eigenvectors, const int64_t (*triplets)[3],
     const int64_t num_triplets, const int64_t (*bz_grid_addresses)[3],
     const int64_t D_diag[3], const int64_t Q[3][3], const double *fc3,
-    const int64_t is_compact_fc3, const double (*svecs)[3],
-    const int64_t multi_dims[2], const int64_t (*multiplicity)[2],
-    const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
-    const int64_t *band_indices, const int64_t symmetrize_fc3_q,
-    const int64_t make_r0_average, const char *all_shortest,
-    const double cutoff_frequency, const int64_t openmp_per_triplets) {
+    const char *fc3_nonzero_indices, const int64_t is_compact_fc3,
+    const double (*svecs)[3], const int64_t multi_dims[2],
+    const int64_t (*multiplicity)[2], const double *masses,
+    const int64_t *p2s_map, const int64_t *s2p_map, const int64_t *band_indices,
+    const int64_t symmetrize_fc3_q, const int64_t make_r0_average,
+    const char *all_shortest, const double cutoff_frequency,
+    const int64_t openmp_per_triplets) {
     RecgridConstBZGrid *bzgrid;
     AtomTriplets *atom_triplets;
     int64_t i, j;
@@ -102,6 +103,7 @@ int64_t ph3py_get_interaction(
     atom_triplets->s2p_map = s2p_map;
     atom_triplets->make_r0_average = make_r0_average;
     atom_triplets->all_shortest = all_shortest;
+    atom_triplets->nonzero_indices = fc3_nonzero_indices;
 
     itr_get_interaction(fc3_normal_squared, g_zero, frequencies,
                         (lapack_complex_double *)eigenvectors, triplets,
@@ -127,7 +129,8 @@ int64_t ph3py_get_pp_collision(
     const int64_t (*bz_grid_addresses)[3], /* thm */
     const int64_t *bz_map,                 /* thm */
     const int64_t bz_grid_type, const int64_t D_diag[3], const int64_t Q[3][3],
-    const double *fc3, const int64_t is_compact_fc3, const double (*svecs)[3],
+    const double *fc3, const char *fc3_nonzero_indices,
+    const int64_t is_compact_fc3, const double (*svecs)[3],
     const int64_t multi_dims[2], const int64_t (*multiplicity)[2],
     const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
     const Larray *band_indices, const Darray *temperatures_THz,
@@ -169,6 +172,7 @@ int64_t ph3py_get_pp_collision(
     atom_triplets->s2p_map = s2p_map;
     atom_triplets->make_r0_average = make_r0_average;
     atom_triplets->all_shortest = all_shortest;
+    atom_triplets->nonzero_indices = fc3_nonzero_indices;
 
     ppc_get_pp_collision(imag_self_energy, relative_grid_address, frequencies,
                          (lapack_complex_double *)eigenvectors, triplets,
@@ -192,13 +196,14 @@ int64_t ph3py_get_pp_collision_with_sigma(
     const int64_t (*triplets)[3], const int64_t num_triplets,
     const int64_t *triplet_weights, const int64_t (*bz_grid_addresses)[3],
     const int64_t D_diag[3], const int64_t Q[3][3], const double *fc3,
-    const int64_t is_compact_fc3, const double (*svecs)[3],
-    const int64_t multi_dims[2], const int64_t (*multiplicity)[2],
-    const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
-    const Larray *band_indices, const Darray *temperatures_THz,
-    const int64_t is_NU, const int64_t symmetrize_fc3_q,
-    const int64_t make_r0_average, const char *all_shortest,
-    const double cutoff_frequency, const int64_t openmp_per_triplets) {
+    const char *fc3_nonzero_indices, const int64_t is_compact_fc3,
+    const double (*svecs)[3], const int64_t multi_dims[2],
+    const int64_t (*multiplicity)[2], const double *masses,
+    const int64_t *p2s_map, const int64_t *s2p_map, const Larray *band_indices,
+    const Darray *temperatures_THz, const int64_t is_NU,
+    const int64_t symmetrize_fc3_q, const int64_t make_r0_average,
+    const char *all_shortest, const double cutoff_frequency,
+    const int64_t openmp_per_triplets) {
     RecgridConstBZGrid *bzgrid;
     AtomTriplets *atom_triplets;
     int64_t i, j;
@@ -232,6 +237,7 @@ int64_t ph3py_get_pp_collision_with_sigma(
     atom_triplets->s2p_map = s2p_map;
     atom_triplets->make_r0_average = make_r0_average;
     atom_triplets->all_shortest = all_shortest;
+    atom_triplets->nonzero_indices = fc3_nonzero_indices;
 
     ppc_get_pp_collision_with_sigma(
         imag_self_energy, sigma, sigma_cutoff, frequencies,
diff --git a/c/phono3py.h b/c/phono3py.h
index e9cebd99..e5751d79 100644
--- a/c/phono3py.h
+++ b/c/phono3py.h
@@ -53,12 +53,13 @@ int64_t ph3py_get_interaction(
     const _lapack_complex_double *eigenvectors, const int64_t (*triplets)[3],
     const int64_t num_triplets, const int64_t (*bz_grid_addresses)[3],
     const int64_t D_diag[3], const int64_t Q[3][3], const double *fc3,
-    const int64_t is_compact_fc3, const double (*svecs)[3],
-    const int64_t multi_dims[2], const int64_t (*multi)[2],
-    const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
-    const int64_t *band_indices, const int64_t symmetrize_fc3_q,
-    const int64_t make_r0_average, const char *all_shortest,
-    const double cutoff_frequency, const int64_t openmp_per_triplets);
+    const char *fc3_nonzero_indices, const int64_t is_compact_fc3,
+    const double (*svecs)[3], const int64_t multi_dims[2],
+    const int64_t (*multi)[2], const double *masses, const int64_t *p2s_map,
+    const int64_t *s2p_map, const int64_t *band_indices,
+    const int64_t symmetrize_fc3_q, const int64_t make_r0_average,
+    const char *all_shortest, const double cutoff_frequency,
+    const int64_t openmp_per_triplets);
 int64_t ph3py_get_pp_collision(
     double *imag_self_energy,
     const int64_t relative_grid_address[24][4][3], /* thm */
@@ -68,7 +69,8 @@ int64_t ph3py_get_pp_collision(
     const int64_t (*bz_grid_addresses)[3], /* thm */
     const int64_t *bz_map,                 /* thm */
     const int64_t bz_grid_type, const int64_t D_diag[3], const int64_t Q[3][3],
-    const double *fc3, const int64_t is_compact_fc3, const double (*svecs)[3],
+    const double *fc3, const char *fc3_nonzero_indices,
+    const int64_t is_compact_fc3, const double (*svecs)[3],
     const int64_t multi_dims[2], const int64_t (*multi)[2],
     const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
     const Larray *band_indices, const Darray *temperatures_THz,
@@ -81,13 +83,14 @@ int64_t ph3py_get_pp_collision_with_sigma(
     const int64_t (*triplets)[3], const int64_t num_triplets,
     const int64_t *triplet_weights, const int64_t (*bz_grid_addresses)[3],
     const int64_t D_diag[3], const int64_t Q[3][3], const double *fc3,
-    const int64_t is_compact_fc3, const double (*svecs)[3],
-    const int64_t multi_dims[2], const int64_t (*multi)[2],
-    const double *masses, const int64_t *p2s_map, const int64_t *s2p_map,
-    const Larray *band_indices, const Darray *temperatures_THz,
-    const int64_t is_NU, const int64_t symmetrize_fc3_q,
-    const int64_t make_r0_average, const char *all_shortest,
-    const double cutoff_frequency, const int64_t openmp_per_triplets);
+    const char *fc3_nonzero_indices, const int64_t is_compact_fc3,
+    const double (*svecs)[3], const int64_t multi_dims[2],
+    const int64_t (*multi)[2], const double *masses, const int64_t *p2s_map,
+    const int64_t *s2p_map, const Larray *band_indices,
+    const Darray *temperatures_THz, const int64_t is_NU,
+    const int64_t symmetrize_fc3_q, const int64_t make_r0_average,
+    const char *all_shortest, const double cutoff_frequency,
+    const int64_t openmp_per_triplets);
 void ph3py_get_imag_self_energy_at_bands_with_g(
     double *imag_self_energy, const Darray *fc3_normal_squared,
     const double *frequencies, const int64_t (*triplets)[3],
diff --git a/c/real_to_reciprocal.c b/c/real_to_reciprocal.c
index cb125528..a5cc46d9 100644
--- a/c/real_to_reciprocal.c
+++ b/c/real_to_reciprocal.c
@@ -277,7 +277,7 @@ static void real_to_reciprocal_elements(
     const int64_t pi0, const int64_t pi1, const int64_t pi2,
     const int64_t leg_index) {
     int64_t i, j, k, l;
-    int64_t num_satom, adrs_shift;
+    int64_t num_satom, adrs_shift_atoms, adrs_shift_fc3;
     lapack_complex_double phase_factor;
     double fc3_rec_real[27], fc3_rec_imag[27];
 
@@ -309,8 +309,11 @@ static void real_to_reciprocal_elements(
                     continue;
                 }
             }
-            adrs_shift =
-                i * 27 * num_satom * num_satom + j * 27 * num_satom + k * 27;
+            adrs_shift_atoms = i * num_satom * num_satom + j * num_satom + k;
+            if (!atom_triplets->nonzero_indices[adrs_shift_atoms]) {
+                continue;
+            }
+            adrs_shift_fc3 = adrs_shift_atoms * 27;
             phase_factor =
                 phonoc_complex_prod(phase_factor1[j], phase_factor2[k]);
 
@@ -320,19 +323,19 @@ static void real_to_reciprocal_elements(
                 for (l = 0; l < 27; l++) {
                     fc3_rec_real[l] +=
                         lapack_complex_double_real(phase_factor) *
-                        fc3[adrs_shift + l] * 3;
+                        fc3[adrs_shift_fc3 + l] * 3;
                     fc3_rec_imag[l] +=
                         lapack_complex_double_imag(phase_factor) *
-                        fc3[adrs_shift + l] * 3;
+                        fc3[adrs_shift_fc3 + l] * 3;
                 }
             } else {
                 for (l = 0; l < 27; l++) {
                     fc3_rec_real[l] +=
                         lapack_complex_double_real(phase_factor) *
-                        fc3[adrs_shift + l];
+                        fc3[adrs_shift_fc3 + l];
                     fc3_rec_imag[l] +=
                         lapack_complex_double_imag(phase_factor) *
-                        fc3[adrs_shift + l];
+                        fc3[adrs_shift_fc3 + l];
                 }
             }
         }
diff --git a/c/real_to_reciprocal.h b/c/real_to_reciprocal.h
index 5919965f..a03cc035 100644
--- a/c/real_to_reciprocal.h
+++ b/c/real_to_reciprocal.h
@@ -48,6 +48,7 @@ typedef struct {
     const int64_t *s2p_map;
     int64_t make_r0_average;
     const char *all_shortest;
+    const char *nonzero_indices;  // for compact fc3
 } AtomTriplets;
 
 void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py
index 3fabc16b..985f3861 100644
--- a/phono3py/api_phono3py.py
+++ b/phono3py/api_phono3py.py
@@ -1138,7 +1138,7 @@ class Phono3py:
             `None` will choose one of them automatically.
 
         """
-        if self.mesh_numbers is None:
+        if self._bz_grid is None:
             msg = "Phono3py.mesh_numbers of instance has to be set."
             raise RuntimeError(msg)
 
@@ -1151,6 +1151,7 @@ class Phono3py:
             self._bz_grid,
             self._primitive_symmetry,
             fc3=self._fc3,
+            fc3_nonzero_indices=self._fc3_nonzero_indices,
             band_indices=self._band_indices_flatten,
             constant_averaged_interaction=constant_averaged_interaction,
             frequency_factor_to_THz=self._frequency_factor_to_THz,
diff --git a/phono3py/conductivity/rta_base.py b/phono3py/conductivity/rta_base.py
index d065457f..2ac45d70 100644
--- a/phono3py/conductivity/rta_base.py
+++ b/phono3py/conductivity/rta_base.py
@@ -328,7 +328,6 @@ class ConductivityRTABase(ConductivityBase):
             s2p,
             masses,
         ) = self._pp.get_primitive_and_supercell_correspondence()
-        fc3 = self._pp.fc3
         triplets_at_q, weights_at_q, _, _ = self._pp.get_triplets_at_q()
 
         if None in self._sigmas:
@@ -384,7 +383,8 @@ class ConductivityRTABase(ConductivityBase):
                     self._pp.bz_grid.store_dense_gp_map * 1 + 1,
                     self._pp.bz_grid.D_diag,
                     self._pp.bz_grid.Q,
-                    fc3,
+                    self._pp.fc3,
+                    self._pp.fc3_nonzero_indices,
                     svecs,
                     multi,
                     masses,
@@ -415,7 +415,8 @@ class ConductivityRTABase(ConductivityBase):
                     self._pp.bz_grid.addresses,
                     self._pp.bz_grid.D_diag,
                     self._pp.bz_grid.Q,
-                    fc3,
+                    self._pp.fc3,
+                    self._pp.fc3_nonzero_indices,
                     svecs,
                     multi,
                     masses,
diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py
index 37b096fc..0edcd2ef 100644
--- a/phono3py/cui/load.py
+++ b/phono3py/cui/load.py
@@ -509,7 +509,7 @@ def _load_fc3(
         assert "fc3_nonzero_indices" in fc3
         ph3py.fc3_nonzero_indices = fc3["fc3_nonzero_indices"]
         if log_level:
-            print(f'fc3 and nonzero indices were read from "{_fc3_filename}".')
+            print(f'fc3 and fc3 nonzero indices were read from "{_fc3_filename}".')
     else:
         _check_fc3_shape(ph3py, fc3, filename=_fc3_filename)
         ph3py.fc3 = fc3
diff --git a/phono3py/phonon3/interaction.py b/phono3py/phonon3/interaction.py
index 76c3affa..4cd27b08 100644
--- a/phono3py/phonon3/interaction.py
+++ b/phono3py/phonon3/interaction.py
@@ -36,9 +36,10 @@
 from __future__ import annotations
 
 from collections.abc import Sequence
-from typing import Literal, Optional
+from typing import Literal
 
 import numpy as np
+from numpy.typing import NDArray
 from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, get_dynamical_matrix
 from phonopy.physical_units import get_physical_units
 from phonopy.structure.cells import Primitive, compute_all_sg_permutations
@@ -93,8 +94,9 @@ class Interaction:
         primitive: Primitive,
         bz_grid: BZGrid,
         primitive_symmetry: Symmetry,
-        fc3: np.ndarray | None = None,
-        band_indices: np.ndarray | Sequence | None = None,
+        fc3: NDArray | None = None,
+        fc3_nonzero_indices: NDArray | None = None,
+        band_indices: NDArray | Sequence | None = None,
         constant_averaged_interaction: float | None = None,
         frequency_factor_to_THz: float | None = None,
         frequency_scale_factor: float | None = None,
@@ -120,7 +122,7 @@ class Interaction:
         self._frequency_scale_factor = frequency_scale_factor
 
         if fc3 is not None:
-            self._set_fc3(fc3)
+            self._set_fc3(fc3, fc3_nonzero_indices=fc3_nonzero_indices)
 
         # Unit to eV^2
         if unit_conversion is None:
@@ -180,9 +182,7 @@ class Interaction:
         )
         self._get_all_shortest()
 
-    def run(
-        self, lang: Literal["C", "Python"] = "C", g_zero: Optional[np.ndarray] = None
-    ):
+    def run(self, lang: Literal["C", "Python"] = "C", g_zero: NDArray | None = None):
         """Run ph-ph interaction calculation."""
         if self._phonon_all_done:
             self.run_phonon_solver()
@@ -209,7 +209,7 @@ class Interaction:
             )
 
     @property
-    def interaction_strength(self) -> np.ndarray | None:
+    def interaction_strength(self) -> NDArray | None:
         """Return ph-ph interaction strength.
 
         Returns
@@ -222,7 +222,7 @@ class Interaction:
         return self._interaction_strength
 
     @property
-    def mesh_numbers(self) -> np.ndarray:
+    def mesh_numbers(self) -> NDArray:
         """Return mesh numbers.
 
         Returns
@@ -239,10 +239,15 @@ class Interaction:
         return self._is_mesh_symmetry
 
     @property
-    def fc3(self) -> np.ndarray:
+    def fc3(self) -> NDArray:
         """Return fc3."""
         return self._fc3
 
+    @property
+    def fc3_nonzero_indices(self) -> NDArray:
+        """Return fc3_nonzero_indices."""
+        return self._fc3_nonzero_indices
+
     @property
     def dynamical_matrix(self) -> DynamicalMatrix | None:
         """Return DynamicalMatrix class instance."""
@@ -273,10 +278,10 @@ class Interaction:
     def get_triplets_at_q(
         self,
     ) -> tuple[
-        np.ndarray | None,
-        np.ndarray | None,
-        np.ndarray | None,
-        np.ndarray | None,
+        NDArray | None,
+        NDArray | None,
+        NDArray | None,
+        NDArray | None,
     ]:
         """Return grid point triplets information.
 
@@ -299,7 +304,7 @@ class Interaction:
         return self._bz_grid
 
     @property
-    def band_indices(self) -> np.ndarray:
+    def band_indices(self) -> NDArray:
         """Return band indices.
 
         Returns
@@ -316,7 +321,7 @@ class Interaction:
         return self._nac_params
 
     @property
-    def nac_q_direction(self) -> np.ndarray | None:
+    def nac_q_direction(self) -> NDArray | None:
         """Return q-direction used for NAC at q->0.
 
         Direction of q-vector watching from Gamma point used for
@@ -337,7 +342,7 @@ class Interaction:
             self._nac_q_direction = np.array(nac_q_direction, copy=True, dtype="double")
 
     @property
-    def zero_value_positions(self) -> np.ndarray | None:
+    def zero_value_positions(self) -> NDArray | None:
         """Return zero ph-ph interaction elements information.
 
         Returns
@@ -349,7 +354,7 @@ class Interaction:
 
     def get_phonons(
         self,
-    ) -> tuple[np.ndarray | None, np.ndarray | None, np.ndarray | None]:
+    ) -> tuple[NDArray | None, NDArray | None, NDArray | None]:
         """Return phonons on grid.
 
         Returns
@@ -407,7 +412,7 @@ class Interaction:
         return self._make_r0_average
 
     @property
-    def all_shortest(self) -> np.ndarray:
+    def all_shortest(self) -> NDArray:
         """Return boolean of make_r0_average.
 
         This flag is used to activate averaging of fc3 transformation
@@ -419,7 +424,7 @@ class Interaction:
         return self._all_shortest
 
     @property
-    def averaged_interaction(self) -> np.ndarray:
+    def averaged_interaction(self) -> NDArray:
         """Return sum over phonon triplets of interaction strength.
 
         See Eq.(21) of PRB 91, 094306 (2015)
@@ -439,7 +444,7 @@ class Interaction:
 
     def get_primitive_and_supercell_correspondence(
         self,
-    ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    ) -> tuple[NDArray, NDArray, NDArray, NDArray, NDArray]:
         """Return atomic pair information."""
         return (self._svecs, self._multi, self._p2s, self._s2p, self._masses)
 
@@ -449,7 +454,7 @@ class Interaction:
         return self._unit_conversion
 
     @property
-    def constant_averaged_interaction(self) -> Optional[float]:
+    def constant_averaged_interaction(self) -> float | None:
         """Return constant averaged interaction."""
         return self._constant_averaged_interaction
 
@@ -817,7 +822,7 @@ class Interaction:
         self._interaction_strength = None
         self._g_zero = None
 
-    def _set_fc3(self, fc3):
+    def _set_fc3(self, fc3: NDArray, fc3_nonzero_indices: NDArray | None = None):
         if (
             isinstance(fc3, np.ndarray)
             and fc3.dtype == np.dtype("double")
@@ -834,7 +839,24 @@ class Interaction:
                 fc3 * self._frequency_scale_factor**2, dtype="double", order="C"
             )
 
-    def _get_band_indices(self, band_indices) -> np.ndarray:
+        if fc3_nonzero_indices is None:
+            self._fc3_nonzero_indices = np.ones(
+                self._fc3.shape[:3], dtype="byte", order="C"
+            )
+        elif (
+            isinstance(fc3_nonzero_indices, np.ndarray)
+            and fc3_nonzero_indices.dtype == np.dtype("byte")
+            and fc3_nonzero_indices.flags.aligned
+            and fc3_nonzero_indices.flags.owndata
+            and fc3_nonzero_indices.flags.c_contiguous
+        ):
+            self._fc3_nonzero_indices = fc3_nonzero_indices
+        else:
+            self._fc3_nonzero_indices = np.array(
+                fc3_nonzero_indices, dtype="byte", order="C"
+            )
+
+    def _get_band_indices(self, band_indices) -> NDArray:
         num_band = len(self._primitive) * 3
         if band_indices is None:
             return np.arange(num_band, dtype="int64")
@@ -877,6 +899,7 @@ class Interaction:
             self._bz_grid.D_diag,
             self._bz_grid.Q,
             self._fc3,
+            self._fc3_nonzero_indices,
             self._svecs,
             self._multi,
             self._masses,
diff --git a/test/conductivity/test_kappa_RTA.py b/test/conductivity/test_kappa_RTA.py
index ea82c728..0932972d 100644
--- a/test/conductivity/test_kappa_RTA.py
+++ b/test/conductivity/test_kappa_RTA.py
@@ -276,6 +276,18 @@ def test_kappa_RTA_nacl(nacl_pbe: Phono3py):
     np.testing.assert_allclose(ref_kappa_RTA, kappa, atol=0.5)
 
 
+def test_kappa_RTA_mgo(
+    mgo_222rd_444rd_symfc: Phono3py, mgo_222rd_444rd_symfc_compact_fc: Phono3py
+):
+    """Test RTA by MgO cutoff 4."""
+    for ph3 in (mgo_222rd_444rd_symfc, mgo_222rd_444rd_symfc_compact_fc):
+        ref_kappa_RTA = [63.75, 63.75, 63.75, 0, 0, 0]
+        kappa = _get_kappa(ph3, [11, 11, 11]).ravel()
+        np.testing.assert_allclose(ref_kappa_RTA, kappa, atol=0.5)
+        kappa = _get_kappa(ph3, [11, 11, 11], is_full_pp=True).ravel()
+        np.testing.assert_allclose(ref_kappa_RTA, kappa, atol=0.5)
+
+
 def test_kappa_RTA_nacl_with_sigma(nacl_pbe: Phono3py):
     """Test RTA with smearing method by NaCl."""
     if nacl_pbe._make_r0_average:
diff --git a/test/conftest.py b/test/conftest.py
index 65917a7f..fa04fc82 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -566,6 +566,48 @@ def mgo_222rd_444rd() -> Phono3py:
     return phono3py.load(yaml_filename, produce_fc=False, log_level=1)
 
 
+@pytest.fixture(scope="session")
+def mgo_222rd_444rd_symfc() -> Phono3py:
+    """Return Phono3py instance of MgO-2x2x2-4x4x4 RD-RD.
+
+    * with symmetry
+    * full fc
+    * use symfc if available on test side
+
+    """
+    pytest.importorskip("symfc")
+
+    yaml_filename = cwd / "phono3py_params_MgO-222rd-444rd.yaml.xz"
+    return phono3py.load(
+        yaml_filename,
+        is_compact_fc=False,
+        fc_calculator="symfc",
+        fc_calculator_options="|cutoff = 4",
+        log_level=1,
+    )
+
+
+@pytest.fixture(scope="session")
+def mgo_222rd_444rd_symfc_compact_fc() -> Phono3py:
+    """Return Phono3py instance of MgO-2x2x2-4x4x4 RD-RD.
+
+    * with symmetry
+    * full fc
+    * use symfc if available on test side
+
+    """
+    pytest.importorskip("symfc")
+
+    yaml_filename = cwd / "phono3py_params_MgO-222rd-444rd.yaml.xz"
+    return phono3py.load(
+        yaml_filename,
+        is_compact_fc=True,
+        fc_calculator="symfc",
+        fc_calculator_options="|cutoff = 4",
+        log_level=1,
+    )
+
+
 @pytest.fixture(scope="session")
 def ph_nacl() -> Phonopy:
     """Return Phonopy class instance of NaCl 2x2x2."""
diff --git a/test/phonon3/test_interaction.py b/test/phonon3/test_interaction.py
index 8af0226b..2a94550d 100644
--- a/test/phonon3/test_interaction.py
+++ b/test/phonon3/test_interaction.py
@@ -301,6 +301,7 @@ def test_get_all_shortest(aln_lda: Phono3py):
     """Test Interaction._get_all_shortest."""
     ph3 = aln_lda
     ph3.mesh_numbers = 30
+    assert ph3.grid is not None
     itr = Interaction(
         ph3.primitive,
         ph3.grid,
@@ -339,6 +340,7 @@ def _get_irt(
     make_r0_average: bool = False,
 ):
     ph3.mesh_numbers = mesh
+    assert ph3.grid is not None
     itr = Interaction(
         ph3.primitive,
         ph3.grid,