diff --git a/.github/workflows/phono3py-pytest-conda-mkl-phphmtblas.yml b/.github/workflows/phono3py-pytest-conda-mkl-phphmtblas.yml
index a262cdf6..49e60dd4 100644
--- a/.github/workflows/phono3py-pytest-conda-mkl-phphmtblas.yml
+++ b/.github/workflows/phono3py-pytest-conda-mkl-phphmtblas.yml
@@ -24,7 +24,7 @@ jobs:
       run: |
         conda activate test
         conda install --yes python=${{ matrix.python-version }}
-        conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest spglib alm cmake c-compiler cxx-compiler
+        conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest spglib alm cmake c-compiler cxx-compiler pypolymlp
     - name: Install symfc develop branch
       run: |
         conda activate test
diff --git a/.github/workflows/phono3py-pytest-conda-mkl-v2.yml b/.github/workflows/phono3py-pytest-conda-mkl-v2.yml
index a6d5c8d5..c1e95df6 100644
--- a/.github/workflows/phono3py-pytest-conda-mkl-v2.yml
+++ b/.github/workflows/phono3py-pytest-conda-mkl-v2.yml
@@ -24,7 +24,7 @@ jobs:
       run: |
         conda activate test
         conda install --yes python=${{ matrix.python-version }}
-        conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest spglib alm cmake c-compiler cxx-compiler
+        conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest spglib alm cmake c-compiler cxx-compiler pypolymlp
     - name: Install symfc develop branch
       run: |
         conda activate test
diff --git a/.github/workflows/phono3py-pytest-conda-mkl.yml b/.github/workflows/phono3py-pytest-conda-mkl.yml
index f42b8264..c1ae62de 100644
--- a/.github/workflows/phono3py-pytest-conda-mkl.yml
+++ b/.github/workflows/phono3py-pytest-conda-mkl.yml
@@ -24,7 +24,7 @@ jobs:
       run: |
         conda activate test
         conda install --yes python=${{ matrix.python-version }}
-        conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest codecov pytest-cov spglib alm cmake c-compiler cxx-compiler
+        conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest codecov pytest-cov spglib alm cmake c-compiler cxx-compiler pypolymlp
     - name: Install symfc develop branch
       run: |
         conda activate test
diff --git a/.github/workflows/phono3py-pytest-conda-numpy2.yml b/.github/workflows/phono3py-pytest-conda-numpy2.yml
index 85aee8be..c2e05745 100644
--- a/.github/workflows/phono3py-pytest-conda-numpy2.yml
+++ b/.github/workflows/phono3py-pytest-conda-numpy2.yml
@@ -24,7 +24,7 @@ jobs:
       run: |
         conda activate test
         conda install --yes python=${{ matrix.python-version }}
-        conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py "numpy=2" scipy pytest cmake c-compiler cxx-compiler
+        conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py "numpy=2" scipy pytest cmake c-compiler cxx-compiler pypolymlp
     - name: Install spglib develop branch
       run: |
         conda activate test
diff --git a/.github/workflows/phono3py-pytest-conda-phphmtblas.yml b/.github/workflows/phono3py-pytest-conda-phphmtblas.yml
index 5771b33b..cd4d625e 100644
--- a/.github/workflows/phono3py-pytest-conda-phphmtblas.yml
+++ b/.github/workflows/phono3py-pytest-conda-phphmtblas.yml
@@ -24,7 +24,7 @@ jobs:
       run: |
         conda activate test
         conda install --yes python=${{ matrix.python-version }}
-        conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py scipy pytest spglib alm cmake c-compiler cxx-compiler
+        conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py scipy pytest spglib alm cmake c-compiler cxx-compiler pypolymlp
     - name: Install symfc develop branch
       run: |
         conda activate test
diff --git a/.github/workflows/phono3py-pytest-conda.yml b/.github/workflows/phono3py-pytest-conda.yml
index 1df4515a..66bb1385 100644
--- a/.github/workflows/phono3py-pytest-conda.yml
+++ b/.github/workflows/phono3py-pytest-conda.yml
@@ -20,11 +20,25 @@ jobs:
     - uses: conda-incubator/setup-miniconda@v3
       with:
         miniforge-version: latest
-    - name: Install dependent packages
+    - name: Install dependent packages for python > 3.9
+      if: ${{ matrix.python-version != 3.9 }}
+      run: |
+        conda activate test
+        conda install --yes python=${{ matrix.python-version }}
+        conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py "numpy>=2.1" scipy pytest spglib alm cmake c-compiler cxx-compiler pypolymlp
+    - name: Install dependent packages for python == 3.9
+      if: ${{ matrix.python-version == 3.9 }}
       run: |
         conda activate test
         conda install --yes python=${{ matrix.python-version }}
         conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py scipy pytest spglib alm cmake c-compiler cxx-compiler
+    - name: Install symfc develop branch
+      run: |
+        conda activate test
+        git clone --depth 1 https://github.com/symfc/symfc.git
+        cd symfc
+        pip install -e . -vvv
+        cd ..
     - name: Install phonopy develop branch
       run: |
         conda activate test
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 94ea09a8..6d3c4221 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -2,7 +2,7 @@
 # See https://pre-commit.com/hooks.html for more hooks
 repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.6.0
+  rev: v5.0.0
   hooks:
   - id: trailing-whitespace
   - id: end-of-file-fixer
@@ -11,7 +11,7 @@ repos:
     exclude: ^example/AlN-LDA/
 
 - repo: https://github.com/astral-sh/ruff-pre-commit
-  rev: v0.6.7
+  rev: v0.7.0
   hooks:
   - id: ruff
     args: [ "--fix", "--show-fixes" ]
diff --git a/doc/changelog.md b/doc/changelog.md
index 71a671f2..4faa2c1f 100644
--- a/doc/changelog.md
+++ b/doc/changelog.md
@@ -2,6 +2,10 @@
 
 # Change Log
 
+## Nov-3-2024: Version 3.6.0
+
+- Maintenance release.
+
 ## Sep-24-2024: Version 3.5.2
 
 - Fix a memory leak.
diff --git a/doc/command-options.md b/doc/command-options.md
index b4957714..33565035 100644
--- a/doc/command-options.md
+++ b/doc/command-options.md
@@ -626,9 +626,9 @@ $$
  \bigl|\Phi_{-\lambda\lambda'\lambda''}\bigl|^2
  \left\{(n_{\lambda'}+ n_{\lambda''}+1)
   \delta(\omega-\omega_{\lambda'}-\omega_{\lambda''}) \right.
-  * (n_{\lambda'}-n_{\lambda''})
+  + (n_{\lambda'}-n_{\lambda''})
  \left[\delta(\omega+\omega_{\lambda'}-\omega_{\lambda''})
-* \left. \delta(\omega-\omega_{\lambda'}+\omega_{\lambda''})
+- \left. \delta(\omega-\omega_{\lambda'}+\omega_{\lambda''})
 \right]\right\}.
 $$
 
@@ -858,7 +858,7 @@ $$
  \bigl|\Phi_{-\lambda\lambda_1\lambda_2}\bigl|^2
  (n_{\lambda_1}-n_{\lambda_2})
  \left[\delta(\omega+\omega_{\lambda_1}-\omega_{\lambda_2})
- *  \delta(\omega-\omega_{\lambda_1}+\omega_{\lambda_2})
+ - \delta(\omega-\omega_{\lambda_1}+\omega_{\lambda_2})
  \right]
  \end{align*}
 $$
@@ -872,7 +872,7 @@ $$
  \bigl|\Phi_{-\lambda\lambda_1\lambda_2}\bigl|^2
 (n_{\lambda_1}+ n_{\lambda_2}+1)
  \left[ \delta(\omega-\omega_{\lambda_1}-\omega_{\lambda_2})
- * \delta(\omega + \omega_{\lambda_1} + \omega_{\lambda_2}) \right]
+ - \delta(\omega + \omega_{\lambda_1} + \omega_{\lambda_2}) \right]
  \end{align*},
 $$
 
@@ -940,11 +940,11 @@ $$
  \bigl|\Phi_{-\lambda\lambda_1\lambda_2}\bigl|^2 &
  \left\{(n_{\lambda_1}+ n_{\lambda_2}+1)
  \left[ \delta(\omega-\omega_{\lambda_1}-\omega_{\lambda_2})
- * \delta(\omega+\omega_{\lambda_1}+\omega_{\lambda_2}) \right] \right.
+ - \delta(\omega+\omega_{\lambda_1}+\omega_{\lambda_2}) \right] \right.
  \\
  & + (n_{\lambda_1}-n_{\lambda_2})
  \left[\delta(\omega+\omega_{\lambda_1}-\omega_{\lambda_2})
- * \left. \delta(\omega-\omega_{\lambda_1}+\omega_{\lambda_2})
+ - \left. \delta(\omega-\omega_{\lambda_1}+\omega_{\lambda_2})
  \right]\right\},
 \end{align*}
 $$
@@ -959,14 +959,14 @@ $$
  \left\{
  \left[ \frac{(n_{\lambda_1}+ n_{\lambda_2}+1)}{
  (\omega-\omega_{\lambda_1}-\omega_{\lambda_2})_\mathrm{p}}
- * \frac{(n_{\lambda_1}+ n_{\lambda_2}+1)}{
+ - \frac{(n_{\lambda_1}+ n_{\lambda_2}+1)}{
  (\omega+\omega_{\lambda_1}+\omega_{\lambda_2})_\mathrm{p}}
  \right]
  \right. \\
   & + \left[
  \frac{(n_{\lambda_1}-n_{\lambda_2})}{(\omega +
  \omega_{\lambda_1} - \omega_{\lambda_2})_\mathrm{p}}
- * \left. \frac{(n_{\lambda_1}-n_{\lambda_2})}{(\omega -
+ - \left. \frac{(n_{\lambda_1}-n_{\lambda_2})}{(\omega -
  \omega_{\lambda_1} + \omega_{\lambda_2})_\mathrm{p}}
  \right]\right\},
 \end{align*}
@@ -1030,7 +1030,7 @@ A_\lambda(\omega) = \frac{1}{\pi} \frac{4\Omega^2_\lambda
  \Gamma_\lambda(\omega)}
  {\left[\omega^2 - \Omega^2_\lambda -
  2\Omega_\lambda \Delta_\lambda(\omega) \right]^2
- * \left[ 2\Omega_\lambda
+ - \left[ 2\Omega_\lambda
  \Gamma_\lambda(\omega) \right]^2},
 $$
 
diff --git a/doc/conf.py b/doc/conf.py
index f186e32a..2f64b5f2 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -58,9 +58,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = "3.5"
+version = "3.6"
 # The full version, including alpha/beta/rc tags.
-release = "3.5.2"
+release = "3.6.0"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/doc/pypolymlp.md b/doc/pypolymlp.md
index 00403198..eb2e2fd3 100644
--- a/doc/pypolymlp.md
+++ b/doc/pypolymlp.md
@@ -23,12 +23,38 @@ For further details on combining phono3py calculations with pypolymlp, refer to
 An example of its usage can be found in the `example/NaCl-pypolymlp` directory
 in the distribution from GitHub or PyPI.
 
+## Citation of pypolymlp
+
+"Tutorial: Systematic development of polynomial machine learning potentials for elemental and alloy systems", A. Seko, J. Appl. Phys. **133**, 011101 (2023).
+
+```
+@article{pypolymlp,
+    author = {Seko, Atsuto},
+    title = "{"Tutorial: Systematic development of polynomial machine learning potentials for elemental and alloy systems"}",
+    journal = {J. Appl. Phys.},
+    volume = {133},
+    number = {1},
+    pages = {011101},
+    year = {2023},
+    month = {01},
+}
+```
+
 ## Requirements
 
 - [pypolymlp](https://github.com/sekocha/pypolymlp)
+
+  For linux (x86-64), a compiled package of pypolymlp can be installed via
+  conda-forge (recommended). Otherwise, pypolymlp can be installed from
+  source-code.
+
 - [symfc](https://github.com/symfc/symfc)
 
-## Workflow
+  Installed via pip, conda-forge, or source code.
+
+## How to calculate
+
+### Workflow
 
 1. Generate random displacements in supercells. Use {ref}`--rd
    <random_displacements_option>` option.
diff --git a/phono3py/cui/create_force_sets.py b/phono3py/cui/create_force_sets.py
index ed5b654f..bcb672ee 100644
--- a/phono3py/cui/create_force_sets.py
+++ b/phono3py/cui/create_force_sets.py
@@ -41,12 +41,18 @@
 from typing import Optional
 
 import numpy as np
-from phonopy.cui.create_force_sets import check_number_of_force_files
+from phonopy.cui.create_force_sets import (
+    check_agreement_of_supercell_positions,
+    check_agreements_of_displacements,
+    check_number_of_force_files,
+)
 from phonopy.cui.load_helper import get_nac_params
 from phonopy.cui.phonopy_script import file_exists, files_exist, print_error
 from phonopy.file_IO import is_file_phonopy_yaml, parse_FORCE_SETS, write_FORCE_SETS
 from phonopy.interface.calculator import get_calc_dataset
+from phonopy.structure.atoms import PhonopyAtoms
 
+from phono3py.cui.settings import Phono3pySettings
 from phono3py.file_IO import (
     get_length_of_first_line,
     parse_FORCES_FC2,
@@ -60,7 +66,7 @@
 
 
 def create_FORCES_FC3_and_FORCES_FC2(
-    settings,
+    settings: Phono3pySettings,
     cell_filename: Optional[str],
     log_level: int = 0,
 ):
@@ -106,11 +112,16 @@ def create_FORCES_FC3_and_FORCES_FC2(
 
     if settings.create_forces_fc3 or settings.create_forces_fc3_file:
         calc_dataset_fc3 = _get_force_sets_fc3(
-            settings, ph3py_yaml.dataset, disp_filename, interface_mode, log_level
+            settings,
+            ph3py_yaml.dataset,
+            ph3py_yaml.supercell,
+            disp_filename,
+            interface_mode,
+            log_level,
         )
         if not calc_dataset_fc3["forces"]:
             if log_level:
-                print("%s could not be created." % "FORCES_FC3")
+                print('"FORCES_FC3" could not be created.')
                 print_error()
             sys.exit(1)
 
@@ -118,13 +129,14 @@ def create_FORCES_FC3_and_FORCES_FC2(
         calc_dataset_fc2 = _get_force_sets_fc2(
             settings,
             ph3py_yaml.phonon_dataset,
+            ph3py_yaml.phonon_supercell,
             disp_filename,
             interface_mode,
             log_level,
         )
         if not calc_dataset_fc2["forces"]:
             if log_level:
-                print("%s could not be created." % "FORCES_FC2")
+                print('"FORCES_FC2" could not be created.')
                 print_error()
             sys.exit(1)
     else:
@@ -157,7 +169,7 @@ def create_FORCES_FC3_and_FORCES_FC2(
                 filename="FORCES_FC3",
             )
             if log_level:
-                print("%s has been created." % "FORCES_FC3")
+                print('"FORCES_FC3" has been created.')
 
         if settings.create_forces_fc2:
             write_FORCES_FC2(
@@ -166,7 +178,7 @@ def create_FORCES_FC3_and_FORCES_FC2(
                 filename="FORCES_FC2",
             )
             if log_level:
-                print("%s has been created." % "FORCES_FC2")
+                print('"FORCES_FC2" has been created.')
 
 
 def create_FORCES_FC2_from_FORCE_SETS(log_level):
@@ -247,7 +259,12 @@ def create_FORCE_SETS_from_FORCES_FCx(
 
 
 def _get_force_sets_fc2(
-    settings, disp_dataset, disp_filename, interface_mode, log_level
+    settings: Phono3pySettings,
+    disp_dataset: dict,
+    supercell: PhonopyAtoms,
+    disp_filename: str,
+    interface_mode: Optional[str],
+    log_level: int,
 ) -> dict:
     interface_mode = settings.calculator
     if log_level:
@@ -277,9 +294,19 @@ def _get_force_sets_fc2(
         verbose=(log_level > 0),
     )
     force_sets = calc_dataset["forces"]
+    if "points" in calc_dataset:
+        if filename := check_agreements_of_displacements(
+            supercell, disp_dataset, calc_dataset["points"], force_filenames
+        ):
+            raise RuntimeError(
+                f'Displacements don\'t match with atomic positions in "{filename}".'
+            )
 
-    if settings.subtract_forces:
-        force_filename = settings.subtract_forces
+    if settings.subtract_forces or settings.subtract_forces_fc2:
+        if settings.subtract_forces_fc2:
+            force_filename = settings.subtract_forces_fc2
+        else:
+            force_filename = settings.subtract_forces
         file_exists(force_filename, log_level=log_level)
         calc_dataset_zero = get_calc_dataset(
             interface_mode,
@@ -289,12 +316,26 @@ def _get_force_sets_fc2(
             ],
             verbose=(log_level > 0),
         )
+        if "points" in calc_dataset_zero:
+            if check_agreement_of_supercell_positions(
+                supercell, calc_dataset_zero["points"][0]
+            ):
+                raise RuntimeError(
+                    "Supercell doesn't match with atomic positions in "
+                    f'"{force_filename}".'
+                )
         force_set_zero = calc_dataset_zero["forces"][0]
         for fs in force_sets:
             fs -= force_set_zero
 
         if log_level > 0:
-            print("Forces in '{force_filename}' were subtracted from supercell forces.")
+            print(
+                f'Forces in "{force_filename}" were subtracted from supercell forces.'
+            )
+
+    if log_level > 0 and "first_atoms" not in disp_dataset:
+        if len(force_filenames) < num_disps:
+            print("** Number of supercell files is less than displacements. **")
 
     if log_level > 0:
         print("")
@@ -303,7 +344,12 @@ def _get_force_sets_fc2(
 
 
 def _get_force_sets_fc3(
-    settings, disp_dataset, disp_filename, interface_mode, log_level
+    settings: Phono3pySettings,
+    disp_dataset: dict,
+    supercell: PhonopyAtoms,
+    disp_filename: str,
+    interface_mode: Optional[str],
+    log_level: int,
 ) -> dict:
     if log_level:
         print(f'FC3 Displacement dataset was read from "{disp_filename}".')
@@ -331,12 +377,14 @@ def _get_force_sets_fc3(
         file_exists(filename, log_level=log_level)
 
     if log_level > 0:
-        print(f"  Number of displacements: {num_disps}")
         print(f"  Number of supercell files: {len(force_filenames)}")
+        print(f'  Number of displacements in "{disp_filename}": {num_disps}')
 
-    if not check_number_of_force_files(num_disps, force_filenames, disp_filename):
+    if "first_atoms" in disp_dataset and not check_number_of_force_files(
+        num_disps, force_filenames, disp_filename
+    ):  # type-1
         calc_dataset = {"forces": []}
-    else:
+    else:  # type-2
         calc_dataset = get_calc_dataset(
             interface_mode,
             num_atoms,
@@ -344,11 +392,18 @@ def _get_force_sets_fc3(
             verbose=(log_level > 0),
         )
         force_sets = calc_dataset["forces"]
+        if "points" in calc_dataset:
+            if filename := check_agreements_of_displacements(
+                supercell, disp_dataset, calc_dataset["points"], force_filenames
+            ):
+                raise RuntimeError(
+                    f'Displacements don\'t match with atomic positions in "{filename}".'
+                )
 
     if settings.subtract_forces:
         force_filename = settings.subtract_forces
         file_exists(force_filename, log_level=log_level)
-        calc_dataset = get_calc_dataset(
+        calc_dataset_zero = get_calc_dataset(
             interface_mode,
             num_atoms,
             [
@@ -356,7 +411,16 @@ def _get_force_sets_fc3(
             ],
             verbose=(log_level > 0),
         )
-        force_set_zero = calc_dataset["forces"][0]
+        force_set_zero = calc_dataset_zero["forces"][0]
+        if "points" in calc_dataset_zero:
+            if check_agreement_of_supercell_positions(
+                supercell, calc_dataset_zero["points"][0]
+            ):
+                raise RuntimeError(
+                    "Supercell doesn't match with atomic positions in "
+                    f'"{force_filename}".'
+                )
+
         for fs in force_sets:
             fs -= force_set_zero
 
@@ -365,6 +429,10 @@ def _get_force_sets_fc3(
                 f"Forces in '{force_filename}' were subtracted from supercell forces."
             )
 
+    if log_level > 0 and "first_atoms" not in disp_dataset:
+        if len(force_filenames) < num_disps:
+            print("** Number of supercell files is less than displacements. **")
+
     if log_level > 0:
         print("")
 
@@ -377,6 +445,7 @@ def _set_forces_and_nac_params(
     calc_dataset_fc3: dict,
     calc_dataset_fc2: Optional[dict],
 ):
+    """Set forces, energies and nac_params to phono3py_yaml."""
     if "first_atoms" in ph3py_yaml.dataset:
         count = len(ph3py_yaml.dataset["first_atoms"])
         for i, d1 in enumerate(ph3py_yaml.dataset["first_atoms"]):
@@ -400,6 +469,12 @@ def _set_forces_and_nac_params(
         ph3py_yaml.dataset["supercell_energies"] = np.array(
             calc_dataset_fc3["supercell_energies"], dtype="double"
         )
+        if len(ph3py_yaml.dataset["forces"]) != len(
+            ph3py_yaml.dataset["displacements"]
+        ):
+            ph3py_yaml.dataset["displacements"] = ph3py_yaml.dataset["displacements"][
+                : len(ph3py_yaml.dataset["forces"])
+            ]
 
     if settings.create_forces_fc2 and calc_dataset_fc2:
         if "first_atoms" in ph3py_yaml.phonon_dataset:
@@ -416,6 +491,12 @@ def _set_forces_and_nac_params(
             ph3py_yaml.phonon_dataset["supercell_energies"] = np.array(
                 calc_dataset_fc2["supercell_energies"], dtype="double"
             )
+            if len(ph3py_yaml.phonon_dataset["forces"]) != len(
+                ph3py_yaml.phonon_dataset["displacements"]
+            ):
+                ph3py_yaml.phonon_dataset["displacements"] = ph3py_yaml.phonon_dataset[
+                    "displacements"
+                ][: len(ph3py_yaml.phonon_dataset["forces"])]
 
     if ph3py_yaml.nac_params is None:
         nac_params = get_nac_params(primitive=ph3py_yaml.primitive)
diff --git a/phono3py/cui/phono3py_argparse.py b/phono3py/cui/phono3py_argparse.py
index 32c392f4..8ac6400f 100644
--- a/phono3py/cui/phono3py_argparse.py
+++ b/phono3py/cui/phono3py_argparse.py
@@ -159,6 +159,14 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
         default=None,
         help="Subtract recidual forces from supercell forces",
     )
+    parser.add_argument(
+        "--cfz-fc2",
+        "--subtract-forces-fc2",
+        metavar="FILE",
+        dest="subtract_forces_fc2",
+        default=None,
+        help="Subtract recidual forces from supercell forces for fc2",
+    )
     parser.add_argument(
         "--cfc",
         "--compact-fc",
diff --git a/phono3py/cui/settings.py b/phono3py/cui/settings.py
index 45f71814..096ee8d2 100644
--- a/phono3py/cui/settings.py
+++ b/phono3py/cui/settings.py
@@ -95,6 +95,7 @@ class Phono3pySettings(Settings):
         "solve_collective_phonon": False,
         "emulate_v2": False,
         "subtract_forces": None,
+        "subtract_forces_fc2": None,
         "temperatures": None,
         "use_ave_pp": False,
         "use_grg": False,
@@ -321,6 +322,10 @@ def set_subtract_forces(self, val):
         """Set subtract_forces."""
         self._v["subtract_forces"] = val
 
+    def set_subtract_forces_fc2(self, val):
+        """Set subtract_forces_fc2."""
+        self._v["subtract_forces_fc2"] = val
+
     def set_temperatures(self, val):
         """Set temperatures."""
         self._v["temperatures"] = val
@@ -605,6 +610,10 @@ def _read_options(self):
             if self._args.subtract_forces:
                 self._confs["subtract_forces"] = self._args.subtract_forces
 
+        if "subtract_forces_fc2" in self._args:
+            if self._args.subtract_forces_fc2:
+                self._confs["subtract_forces_fc2"] = self._args.subtract_forces_fc2
+
         if "temperatures" in self._args:
             if self._args.temperatures is not None:
                 self._confs["temperatures"] = " ".join(self._args.temperatures)
@@ -716,6 +725,7 @@ def _parse_conf(self):
                 "create_forces_fc3_file",
                 "output_yaml_filename",
                 "subtract_forces",
+                "subtract_forces_fc2",
             ):
                 self.set_parameter(conf_key, confs[conf_key])
 
@@ -1009,6 +1019,10 @@ def _set_settings(self):
         if "subtract_forces" in params:
             self._settings.set_subtract_forces(params["subtract_forces"])
 
+        # Subtract residual forces to create FORCES_FC2
+        if "subtract_forces_fc2" in params:
+            self._settings.set_subtract_forces_fc2(params["subtract_forces_fc2"])
+
         # Temperatures for scatterings
         if "temperatures" in params:
             self._settings.set_temperatures(params["temperatures"])
diff --git a/phono3py/version.py b/phono3py/version.py
index d4df8bf1..ea0e6101 100644
--- a/phono3py/version.py
+++ b/phono3py/version.py
@@ -34,4 +34,4 @@
 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
 
-__version__ = "3.5.2"
+__version__ = "3.6.0"
diff --git a/pyproject.toml b/pyproject.toml
index af837355..787f70be 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -16,7 +16,7 @@ dependencies = [
     "matplotlib>=2.2.2",
     "h5py>=3.0",
     "spglib>=2.3",
-    "phonopy>=2.28,<2.29",
+    "phonopy>=2.29,<2.30",
 ]
 license = { file = "LICENSE" }