diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 3d9f8110..dd851ddb 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -7,8 +7,8 @@ repos:
   - id: trailing-whitespace
   - id: end-of-file-fixer
   - id: check-yaml
-    exclude: ^conda/
   - id: check-added-large-files
+    exclude: ^example/AlN-LDA/
 
 - repo: https://github.com/astral-sh/ruff-pre-commit
   rev: v0.5.6
diff --git a/doc/auxiliary-tools.md b/doc/auxiliary-tools.md
index c54c9741..52a44915 100644
--- a/doc/auxiliary-tools.md
+++ b/doc/auxiliary-tools.md
@@ -39,7 +39,7 @@ Let's compute lattice thermal conductivity of Si using the `Si-PBEsol` example
 found in the example directory.
 
 ```bash
-% phono3py --mesh="11 11 11" --sym-fc --br
+% phono3py --mesh="11 11 11" --fc-symmetry --br
 ```
 
 Then using the output file, `kappa-m111111.hdf5`, run `phono3py-kaccum` as
diff --git a/doc/changelog.md b/doc/changelog.md
index 6f6212fe..661b84f6 100644
--- a/doc/changelog.md
+++ b/doc/changelog.md
@@ -2,7 +2,12 @@
 
 # Change Log
 
-## Jul-22-2024: Version 3.3.3
+## Aug-8-2024: Version 3.3.4
+
+- Fix of command line user interface.
+- Fix of phono3py yaml parser.
+
+## Aug-6-2024: Version 3.3.3
 
 - Provide functionality by `--cf3` and `--cf2` command options to create force
   constants from displacement-force dataset of random displacements when
diff --git a/doc/command-options.md b/doc/command-options.md
index 8979cc06..635150fb 100644
--- a/doc/command-options.md
+++ b/doc/command-options.md
@@ -2,6 +2,10 @@
 
 # Command options / Setting tags
 
+(use_config_with_option)=
+
+## Use of configuration file
+
 Phono3py is operated with command options or with a configuration file that
 contains setting tags. In this page, the command options are explained. Most of
 command options have their respective setting tags.
@@ -27,13 +31,19 @@ CELL_FILENAME = POSCAR-unitcell
 where the setting tag names are case insensitive. This is run by
 
 ```bash
-% phono3py setting.conf [command options]
+% phono3py setting.conf [OPTIONS]
 ```
 
 or
 
 ```bash
-% phono3py [command options] -- setting.conf
+% phono3py [OPTIONS] -- setting.conf
+```
+
+When using `phono3py-load` (see also {ref}`phono3py_load_command`)
+
+```bash
+% phono3py-load --config setting.conf [OPTIONS]
 ```
 
 ```{contents}
@@ -201,15 +211,32 @@ When using with `--cf2`, `--cf3` has to be specified simultaneously as below,
 
 ### `--dim` (`DIM`)
 
+**`phono3py-load` doesn't have this option.**
+
 Supercell dimension is specified. See the detail at
-<http://phonopy.github.io/phonopy/setting-tags.html#dim>. When a proper
-`phono3py_disp.yaml` exists in the current directory, this is unnecessary to be
-specified.
+<http://phonopy.github.io/phonopy/setting-tags.html#dim>. When
+`phono3py_disp.yaml` is found in the current directory, it is read
+automatically. Since supercell dimension is written in this file, `--dim` is
+unnecessary to specify. For example, just
+
+```bash
+% phono3py --fc-symmetry
+```
+
+or
+
+```bash
+% phono3py --symfc
+```
+
+can be used to calculate force constants.
 
 (dim_fc2_option)=
 
 ### `--dim-fc2` (`DIM_FC2`)
 
+**`phono3py-load` doesn't have this option.**
+
 Supercell dimension for 2nd order force constants (for harmonic phonons) is
 specified. This is optional. When a proper `phono3py_disp.yaml` exists in the
 current directory, this is unnecessary to be specified.
@@ -228,7 +255,7 @@ supercell size and these force calculations have to be done in addition to the
 usual force calculations for 3rd order force constants.
 
 ```bash
-% phono3py -d --dim="2 2 2" --dim-fc2="4 4 4" -c POSCAR-unitcell
+% phono3py -d --dim 2 2 2 --dim-fc2 4 4 4 --pa auto -c POSCAR-unitcell
 ```
 
 After the force calculations, `--cf2` option is used to create `FORCES_FC2`.
@@ -253,10 +280,11 @@ that created in the usual phono3py run without `--dim-fc2` option.
 ### `--pa`, `--primitive-axes` (`PRIMITIVE_AXES`)
 
 Transformation matrix from a non-primitive cell to the primitive cell. See
-phonopy `PRIMITIVE_AXES` tag (`--pa` option) at
-<http://phonopy.github.io/phonopy/setting-tags.html#primitive-axis>. When a proper
-`phono3py_disp.yaml` exists in the current directory, this is unnecessary to be
-specified.
+phonopy `PRIMITIVE_AXES` tag (`--pa` option) at [primitive-axis
+(phonopy)](https://phonopy.github.io/phonopy/setting-tags.html#primitive-axes-or-primitive-axis).
+When `phono3py_disp.yaml` contains this information and `phono3py_disp.yaml` is
+read when running `phono3py` or `phono3py-load` command, this is unnecessary to
+be specified.
 
 ### `--mass` (`MASS`)
 
@@ -273,24 +301,45 @@ web page](https://phonopy.github.io/phonopy/setting-tags.html#magmom).
 (create_displacements_option)=
 ### `-d` (`CREATE_DISPLACEMENTS = .TRUE.`)
 
-Supercell with displacements are created. Using with `--amplitude` option,
-atomic displacement distances are controlled. With this option, files for
-supercells with displacements and `phono3py_disp.yaml` file are created. `--pa`
-should be specified if the input unit cell structure is not a primitive cell,
-e.g., `--pa="F"` if the input unit cell has F-centring.
+**`phono3py-load` doesn't have this option.**
+
+Supercells with displacements and `phono3py_disp.yaml` are created. Using with
+`--amplitude` option, atomic displacement distances are controlled. With this
+option, files for supercells with displacements and `phono3py_disp.yaml` file
+are created.
 
+It is recommended to use this option with `--pa auto` option to store
+information about primitive cell (`primitive_matrix` key) in
+`phono3py_disp.yaml`, e.g.,
+
+```bash
+% phono3py -c POSCAR-unitcell -d --dim 2 2 2 --dim-fc2 4 4 4 --pa auto
+```
 
 (random_displacements_option)=
 ### `--rd` (`RANDOM_DISPLACEMENTS`), `--rd-fc2` (`RANDOM_DISPLACEMENTS_FC2`) and `--random-seed` (`RANDOM_SEED`)
 
+**`phono3py-load` doesn't have this option.**
+
+See also {ref}`random-displacements`.
+
 Random directional displacements are generated for fc3 and fc2 supercells by
 `--rd` and `--rd-fc2`, respectively. `--amplitude` and `--random-seed` options
 may be used together. These are used in the equivalent way to [`--rd` of
 phonopy](https://phonopy.github.io/phonopy/setting-tags.html#random-displacements).
 
+Like `-d` option, it is recommended to specify `--pa auto` together with `--rd`
+and/or `--rd-fc2`,
+
+```bash
+% phono3py -c POSCAR-unitcell --dim 2 2 2 --dim-fc2 4 4 4 --rd 100 --rd-fc2 2 --pa auto
+```
+
 (amplitude_option)=
 ### `--amplitude` (`DISPLACEMENT_DISTANCE`)
 
+**`phono3py-load` doesn't have this option.**
+
 Atomic displacement distance is specified. This value may be increased for the
 weak interaction systems and decreased when the force calculator is numerically
 very accurate.
@@ -303,13 +352,13 @@ The default value depends on calculator. See
 
 Choice of force constants calculator.
 
-```
-% phono3py --fc-calc symfc ...
+```bash
+% phono3py-load --fc-calc symfc ...
 ```
 
 To use different force constants calculators for fc2 and fc3
-```
-% phono3py --fc-calc "symfc|" ...
+```bash
+% phono3py-load --fc-calc "symfc|" ...
 ```
 Those for fc2 and fc3 are seprated by `|` such as `symfc|` . Blank means to
 employ the finite difference method for systematic displacements generated by
@@ -320,8 +369,8 @@ the option `-d`.
 
 Special options for force constants calculators.
 
-```
-% phono3py --fc-calc-opt "cutoff=8" ...
+```bash
+% phono3py-load --fc-calc-opt "cutoff=8" ...
 ```
 
 Similarly to `--fc-calc`, `|` can be used to separated those for fc2 and fc3.
@@ -336,6 +385,9 @@ Similarly to `--fc-calc`, `|` can be used to separated those for fc2 and fc3.
 
 These are shortcuts of `--fc-calc symfc` and `--fc-calc alm`, respectively.
 
+Please be careful that `--symfc` and `--sym-fc` (deprecated) are similar, but
+different.
+
 ## Force constants
 
 (compact_fc_option)=
@@ -352,16 +404,16 @@ data size becomes large. If the input crystal structure has centring
 {ref}`--pa <pa_option>` is necessary to have smallest data size. In this case,
 `--pa` option has to be specified on reading. Otherwise phono3py can recognize
 if `fc2.hdf5` and `fc3.hdf5` are compact or full automatically. When using with
-`--sym-fc`, the calculated results will become slightly different due to
+`--fc-symmetry`, the calculated results will become slightly different due to
 imperfect symmetrization scheme that phono3py employs.
 
 ```bash
-% phono3py --dim="2 2 2" --cfc --pa="F" -c POSCAR-unitcell
+% phono3py-load --compact-fc
 ```
 
 (symmetrization_option)=
 
-### `--sym-fc` (`FC_SYMMETRY = .TRUE.`)
+### `--fc-symmetry` (`FC_SYMMETRY = .TRUE.`)
 
 Second- and third-order force constants are symmetrized. The index exchange of
 real space force constants and translational invariance symmetry are applied in a
@@ -390,15 +442,16 @@ supercell size and the second choice is using `--cutoff-pair` option.
 
 ### `--cutoff-pair` or `--cutoff-pair-distance` (`CUTOFF_PAIR_DISTANCE`)
 
-This option works differently for the `-d` and `--rd` options.
+This option works in two ways.
 
-For `-d`, A cutoff pair-distance in a supercell is used to reduce the number of
-necessary supercells with displacements to obtain third order force constants.
-As the drawback, a certain number of third-order-force-constants elements are
-abandoned or computed with less numerical accuracy. More details are found at
-{ref}`command_cutoff_pair`.
+When using with `-d` options, a cutoff pair-distance in a supercell is used to
+reduce the number of necessary supercells with displacements to obtain third
+order force constants. As the drawback, a certain number of
+third-order-force-constants elements are abandoned or computed with less
+numerical accuracy. More details are found at {ref}`command_cutoff_pair`.
 
-For `--rd`, `--cutoff-pair VAL` is equivalent to `--fc-calc-opt "cutoff=VAL"`.
+When using with an external force constants calculator, `--cutoff-pair VAL` works
+equivalent to `--fc-calc-opt "cutoff=VAL"`.
 
 ### `--alm`
 
@@ -436,10 +489,10 @@ The mapping table between grid points to its indices is obtained by running with
 `--loglevel=2` option.
 
 ```bash
-% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --write-gamma --gp="0 1 2 3 4 5"
+% phono3py-load --mesh 19 19 19 --br --write-gamma --gp 0 1 2 3 4 5
 ```
 
-where `--gp="0 1 2 3 4 5"` can be also written `--gp="0,1,2,3,4,5"`. `--ga`
+where `--gp 0 1 2 3 4 5` can be also written `--gp="0,1,2,3,4,5"`. `--ga`
 option below can be used similarly for the same purpose.
 
 (ga_option)=
@@ -447,10 +500,10 @@ option below can be used similarly for the same purpose.
 ### `--ga` (`GRID_ADDRESSES`)
 
 This is used to specify grid points like `--gp` option but in their addresses
-represented by integer numbers. For example with `--mesh="16 16 16"`, a q-point
-of (0.5, 0.5, 0.5) is given by `--ga="8 8 8"`. The values have to be integers.
+represented by integer numbers. For example with `--mesh 16 16 16`, a q-point
+of (0.5, 0.5, 0.5) is given by `--ga 8 8 8`. The values have to be integers.
 If you want to specify the point on a path,
-`--ga="0 0 0 1 1 1 2 2 2 3 3 3 ..."`, where each three values are recognized as
+`--ga 0 0 0 1 1 1 2 2 2 3 3 3 ...`, where each three values are recognized as
 a grid point. The grid points given by `--ga` option are translated to grid
 point indices as given by `--gp` option, and the values given by `--ga` option
 will not be shown in log files.
@@ -466,7 +519,7 @@ band indices where the values are calculated and summed and averaged over those
 bands.
 
 ```bash
-% phono3py --fc3 --fc2 --dim="2 2 2" --mesh="16 16 16" -c POSCAR-unitcell --nac --gp="34" --bi="4 5, 6"
+% phono3py-load --mesh 16 16 16 --nac --gp 34 --bi "4 5, 6"
 ```
 
 This option may be also useful to distribute the computational demand such like
@@ -489,7 +542,7 @@ calculated divided these integers by sampling mesh numbers for respective
 reciprocal axes.
 
 ```bash
-% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --wgp
+% phono3py-load --mesh 19 19 19 --wgp
 ```
 
 (stp_option)=
@@ -502,7 +555,7 @@ large a calculation is. Only those for specific grid points are shown by using
 with `--gp` or `--ga` option.
 
 ```bash
-% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --stp --gp 20
+% phono3py-load --mesh 19 19 19 --stp --gp 20
 ```
 
 ## Brillouin zone integration
@@ -533,7 +586,7 @@ This is used when we want to test several $\sigma$ values simultaneously.
 The tails of the Gaussian functions that are used to replace delta functions in
 the equation shown at {ref}`--full-pp <full_pp_option>` are cut with this
 option. The value is specified in number of standard deviation.
-`--sigma-cutoff=5` gives the Gaussian functions to be cut at $5\sigma$. Using
+`--sigma-cutoff 5` gives the Gaussian functions to be cut at $5\sigma$. Using
 this option scarifies the numerical accuracy. So the number has to be carefully
 tested. But computation of phonon-phonon interaction strength becomes much
 faster in exchange for it.
@@ -614,27 +667,27 @@ database of the natural abundance data for elements, which refers Laeter _et
 al._, Pure Appl. Chem., **75**, 683 (2003).
 
 ```bash
-% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br --isotope
+% phono3py-load -v --mesh 32 32 20 --br --isotope
 ```
 
 ### `--mass-variances` or `--mv` (`MASS_VARIANCES`)
 
 Mass variance parameters are specified by this option to include phonon-isotope
 scattering effect in the same way as `--isotope` option. For example of GaN,
-this may be set like `--mv="1.97e-4 1.97e-4 0 0"`. The number of elements has to
+this may be set like `--mv 1.97e-4 1.97e-4 0 0`. The number of elements has to
 correspond to the number of atoms in the primitive cell.
 
 Isotope effect to thermal conductivity may be checked first running without
 isotope calculation:
 
 ```bash
-% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br
+% phono3py-load -v --mesh 32 32 20 --br
 ```
 
 Then running with isotope calculation:
 
 ```bash
-% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br --read-gamma --mv="1.97e-4 1.97e-4 0 0"
+% phono3py -v --mesh 32 32 20 --br --read-gamma --mv 1.97e-4 1.97e-4 0 0
 ```
 
 In the result hdf5 file, currently isotope scattering strength is not written
@@ -648,7 +701,7 @@ A most simple phonon boundary scattering treatment is included. $v_g/L$ is just
 used as the scattering rate, where $v_g$ is the group velocity and $L$ is the
 boundary mean free path. The value is given in micrometre. The default value, 1
 metre, is just used to avoid divergence of phonon lifetime and the contribution
-to the thermal conducitivity is considered negligible.
+to the thermal conductivity is considered negligible.
 
 (ave_pp_option)=
 
@@ -691,13 +744,13 @@ option is strongly recommended.
 First, run full conductivity calculation,
 
 ```bash
-% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br
+% phono3py-load -v --mesh 32 32 20 --br
 ```
 
 Then
 
 ```bash
-% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br --read-gamma --ave-pp -o ave_pp
+% phono3py-load -v --mesh 32 32 20 --br --read-gamma --ave-pp -o ave_pp
 ```
 
 ### `--const-ave-pp` (`CONST_AVE_PP = .TRUE.`)
@@ -715,7 +768,7 @@ $\text{eV}^2$.
 See also {ref}`reference papers <ave_pp_reference>`.
 
 ```bash
-% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br --const-ave-pp=1e-10
+% phono3py-load -v --mesh 32 32 20 --br --const-ave-pp 1e-10
 ```
 
 (normal_umklapp_option)=
@@ -739,7 +792,7 @@ The data are stored in `kappa-mxxx(-gx-sx-sdx).hdf5` file and accessed by
 below:
 
 ```bash
-% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="11 11 11" --fc3 --fc2 --br --nu
+% phono3py-load --mesh 11 11 11 --fc3 --fc2 --br --nu
 ...
 % ipython
 
@@ -820,7 +873,7 @@ $$
 Specific temperatures are specified by `--ts`.
 
 ```bash
-% phono3py --fc3 --fc2 --dim="2 2 2" -v --mesh="11 11 11" -c POSCAR-unitcell --br --ts="200 300 400"
+% phono3py-load -v --mesh 11 11 11 -c POSCAR-unitcell --br --ts 200 300 400
 ```
 
 ### `--tmax`, `--tmin`, `--tstep` (`TMAX`, `TMIN`, `TSTEP`)
@@ -830,7 +883,7 @@ See phonopy's document for the same tags at
 <http://phonopy.github.io/phonopy/setting-tags.html#tprop-tmin-tmax-tstep>.
 
 ```bash
-% phono3py --fc3 --fc2 --dim="2 2 2" -v --mesh="11 11 11" -c POSCAR-unitcell --br --tmin=100 --tmax=1000 --tstep=50
+% phono3py-load -v --mesh 11 11 11 --br --tmin 100 --tmax 1000 --tstep 50
 ```
 
 ## Non-analytical term correction
@@ -918,7 +971,7 @@ to `gammas-mxxx-gx(-sx)-tx-bx.dat` in THz (without $2\pi$) with respect to
 samplied frequency points of $\omega$ in THz (without $2\pi$).
 
 ```bash
-% phono3py --fc3 --fc2 --dim="2 2 2" --mesh="16 16 16" -c POSCAR-unitcell --nac --q-direction="1 0 0" --gp=0 --ise --bi="4 5, 6"
+% phono3py-load --mesh 16 16 16 --nac --q-direction 1 0 0 --gp 0 --ise --bi "4 5, 6"
 ```
 
 (rse_option)=
@@ -945,7 +998,7 @@ $\Delta_\lambda(\omega)$ is written to `deltas-mxxx-gx-sx-tx-bx.dat` in THz
 (without $2\pi$).
 
 ```bash
-% phono3py --fc3 --fc2 --dim="2 2 2" --mesh="16 16 16" -c POSCAR-unitcell --nac --q-direction="1 0 0" --gp=0 --rse --sigma="0.1" --bi="4 5, 6"
+% phono3py-load --mesh 16 16 16 --nac --q-direction 1 0 0 --gp 0 --rse --sigma 0.1 --bi "4 5, 6"
 ```
 
 (spectral_function_option)=
@@ -982,7 +1035,7 @@ THz (without $2\pi$) with respect to samplied frequency points of $\omega$ in
 THz (without $2\pi$), and `spectral-mxxx-gx.hdf5`.
 
 ```bash
-% phono3py --fc3 --fc2 --dim="2 2 2" --mesh="16 16 16" -c POSCAR-unitcell --nac --q-direction="1 0 0" --gp=0 --spf
+% phono3py-load --mesh 16 16 16 --nac --q-direction 1 0 0 --gp 0 --spf
 ```
 
 ```{note}
@@ -1026,7 +1079,7 @@ $$
 See also {ref}`reference papers <spectral_function_reference>`.
 
 ```bash
-% phono3py --fc2 --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="16 16 16" --jdos --ga="0 0 0  8 8 8"
+% phono3py-load --mesh 16 16 16 --jdos --ga 0 0 0 8 8 8
 ```
 
 When temperatures are specified, two classes of weighted JDOS are calculated.
@@ -1052,7 +1105,7 @@ $$
 See also {ref}`reference papers <spectral_function_reference>`.
 
 ```bash
-% phono3py --fc2 --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="16 16 16" --jdos --ga="0 0 0  8 8 8" --ts=300
+% phono3py-load --mesh 16 16 16 --jdos --ga 0 0 0 8 8 8 --ts 300
 ```
 
 This is an example of `Si-PBEsol`.
@@ -1081,13 +1134,13 @@ Mode-Gruneisen-parameters are calculated from fc3.
 Mesh sampling mode:
 
 ```bash
-% phono3py --fc3 --fc2 --dim="2 2 2" -v --mesh="16 16 16" -c POSCAR-unitcell --nac --gruneisen
+% phono3py-load -v --mesh 16 16 16 --nac --gruneisen
 ```
 
 Band path mode:
 
 ```bash
-% phono3py --fc3 --fc2 --dim="2 2 2" -v -c POSCAR-unitcell --nac --gruneisen --band="0 0 0  0 0 1/2"
+% phono3py-load -v --nac --gruneisen --band "0 0 0  0 0 1/2"
 ```
 
 ## File I/O
@@ -1199,7 +1252,7 @@ Phonon frequencies, eigenvectors, and grid point addresses are stored in
 may be required depending on calculation setting.
 
 ```bash
-% phono3py --fc2 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --write-phoonon
+% phono3py-load --mesh 11 11 11 --nac --write-phoonon
 ```
 
 Contents of `phonon-mxxx.hdf5` are watched by:
@@ -1248,7 +1301,7 @@ different CPU architectures. {ref}`--pa <pa_option>` and
 {ref}`--nac <nac_option>` may be required depending on calculation setting.
 
 ```bash
-% phono3py --fc2 --fc3 --dim="2 2  2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --read-phoonon --br
+% phono3py-load --mesh 11 11 11 --nac --read-phoonon --br
 ```
 
 (write_read_pp_option)=
@@ -1266,11 +1319,11 @@ calculation, in writing and reading, ph-ph interaction strength has to be stored
 in memory, so there is overhead in memory than usual RTA calculation.
 
 ```bash
-% phono3py --fc2 --fc3 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --write-pp --br --gp=1
+% phono3py-load --mesh 11 11 11 --nac --write-pp --br --gp 1
 ```
 
 ```bash
-% phono3py --fc2 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --read-pp --br --gp=1
+% phono3py-load --mesh 11 11 11 --nac --read-pp --br --gp 1
 ```
 
 (hdf5_compression_option)=
diff --git a/doc/conf.py b/doc/conf.py
index da5cb499..19a19a0b 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -60,7 +60,7 @@
 # The short X.Y version.
 version = "3.3"
 # The full version, including alpha/beta/rc tags.
-release = "3.3.3"
+release = "3.3.4"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/doc/crystal.md b/doc/crystal.md
index 6cdd00c6..cd81a815 100644
--- a/doc/crystal.md
+++ b/doc/crystal.md
@@ -78,10 +78,10 @@ so the -c crystal.o parameter is not needed.
 4) Calculate 3rd and 2nd order force constants in files `fc3.hdf5` and `fc2.hdf5`:
 
    ```bash
-   % phono3py --crystal --dim="2 2 2" --dim-fc2="4 4 4" --sym-fc
+   % phono3py --crystal --dim="2 2 2" --dim-fc2="4 4 4" --fc-symmetry
    ```
 
-   `--sym-fc` is used to symmetrize second- and third-order force constants.
+   `--fc-symmetry` is used to symmetrize second- and third-order force constants.
 
 5) Thermal conductivity calculation:
 
diff --git a/doc/cutoff-pair.md b/doc/cutoff-pair.md
index a07b1332..88ab4299 100644
--- a/doc/cutoff-pair.md
+++ b/doc/cutoff-pair.md
@@ -1,17 +1,14 @@
 (command_cutoff_pair)=
 
 # Force constants calculation with cutoff pair-distance
-
-Since this calculation is a little bit tricky. It may be recommended to try
+```{note}
+Since usage of this calculation is complicated. It is recommended to try
 {ref}`random-displacements` with `--fc-calc-opt "cutoff = VAL"` before trying
 this option.
+```
 
 Here the detail of the command option {ref}`--cutoff-pair <cutoff_pair_option>`
-is explained.
-
-<!-- See also reference {cite}`Mizokami-PRB-2018`. -->
-
-See also {ref}`a reference paper <cutoff_pair_reference>`.
+is explained. See also {ref}`a reference paper <cutoff_pair_reference>`.
 
 ```{contents}
 :depth: 2
@@ -329,7 +326,7 @@ number_of_pairs_in_cutoff: 110
   300.0     119.501    119.501    119.501     -0.000     -0.000      0.000
   300.0     119.483    119.483    119.483     -0.000     -0.000      0.000
   300.0     119.481    119.481    119.481     -0.000     -0.000      0.000
-% for i in {2..10};do cp phono3py_disp.$i.yaml phono3py_disp.yaml; phono3py --sym-fc --mesh="11 11 11" --br|tee std.sym-$i.out;done
+% for i in {2..10};do cp phono3py_disp.$i.yaml phono3py_disp.yaml; phono3py --fc-symmetry --mesh="11 11 11" --br|tee std.sym-$i.out;done
 % for i in {2..10};do egrep '^\s+300' std.sym-$i.out;done
   300.0     124.650    124.650    124.650     -0.000     -0.000      0.000
   300.0     119.765    119.765    119.765     -0.000     -0.000      0.000
diff --git a/doc/index.md b/doc/index.md
index 1e077455..0ec63523 100644
--- a/doc/index.md
+++ b/doc/index.md
@@ -20,28 +20,16 @@ the Wigner transport equation <LTC_options>`
 Papers that may introduce phono3py:
 
 - Theoretical background is summarized in this paper:
-  <http://dx.doi.org/10.1103/PhysRevB.91.094306> (arxiv
-  <http://arxiv.org/abs/1501.00691>).
+  [PRB.91.094306](http://dx.doi.org/10.1103/PhysRevB.91.094306) (arxiv
+  [1501.00691](http://arxiv.org/abs/1501.00691>)).
 - Introduction to phono3py application:
-  <https://journals.jps.jp/doi/10.7566/JPSJ.92.012001> (open access),
+  [JPSJ.92.012001](https://journals.jps.jp/doi/10.7566/JPSJ.92.012001) (open access),
   and phono3py inputs for 103 compounds found in Fig.17
   <https://github.com/atztogo/phonondb/blob/main/mdr/phono3py_103compounds_fd/README.md>
 - Implementation of phono3py:
-  <https://iopscience.iop.org/article/10.1088/1361-648X/acd831> (open access)
+  [JPCM.35.353001](https://iopscience.iop.org/article/10.1088/1361-648X/acd831)
+  (open access)
 
-A set of pre-calculated data for phono3py in
-
-```{image} Si-kaccum.png
-:width: 20%
-```
-
-```{image} Si-kaccum-MFP.png
-:width: 20%
-```
-
-```{image} Si-kdeplot.png
-:width: 22%
-```
 
 ## Documentation
 
diff --git a/doc/phono3py-load.md b/doc/phono3py-load.md
index 98501bb4..c7f16221 100644
--- a/doc/phono3py-load.md
+++ b/doc/phono3py-load.md
@@ -8,31 +8,42 @@ module. The main aim of introducing this command is to provide uniform usage
 over many different force calculators. Once `phono3py_disp.yaml` is created, the
 following operations will be the same using this command.
 
-The following default behaviours are different from that of those of `phono3py`
-command:
+This is used almost in the same way as `phono3py` command, e.g., but there are
+some differences. The following default behaviours are different from that of
+those of `phono3py` command:
 
 1. `phono3py_xxx.yaml` type file is always necessary in either of two ways:
 
    - `phono3py_xxx.yaml` type file is given as the first argument of the
-     command.
-   - `phono3py_xxx.yaml` type file is put in the current directory with one of
-     the default filenames of `phono3py_params.yaml`, `phono3py_disp.yaml`,
-     `phono3py.yaml`. The searching preference order is `phono3py_params.yaml` >
-     `phono3py_disp.yaml` > `phono3py.yaml`.
+     command, e.g.,
+
+     ```bash
+     % phono3py-load phono3py_xxx.yaml --br --ts 300 --mesh 50
+     ```
+
+   - With first argument unspecified, `phono3py_disp.yaml` or `phono3py.yaml`
+     file is read if it is found in the current directory. If both found,
+     `phono3py_disp.yaml` is read. For example, having `phono3py_disp.yaml`
+     under the current directory,
+
+     ```bash
+     % phono3py-load --br --ts 300 --mesh 50
+     ```
 
 2. `-c` option (read crystal structure) does not exist.
 
-3. `-d` option (create displacements) does not exist.
+3. `-d` option (create displacements) does not exist. Please use `phono3py`
+   command.
 
-4. Use of command options is recommended, but phono3py configuration file can be
-   read through `--config` option.
+4. Phono3py configuration file can be read through `--config` option. See
+   {ref}`use_config_with_option`.
 
 5. If parameters for non-analytical term correction (NAC) are found, NAC is
    automatically enabled. This can be disabled by `--nonac` option.
 
 6. When force constants are calculated from displacements and forces dataset,
    force constants are automatically symmetrized. To disable this, `--no-sym-fc`
-   option is used.
+   option can be used.
 
 7. `-o` option works differently from `phono3py` command. This option requires
    one argument of string. The string is used as the output yaml filename that
diff --git a/doc/qe.md b/doc/qe.md
index db34cce7..537e8eaf 100644
--- a/doc/qe.md
+++ b/doc/qe.md
@@ -67,10 +67,10 @@ only limited number of keywords that are shown in the phonopy web site
    `fc3.hdf5` and `fc2.hdf5` files are created by:
 
    ```bash
-   % phono3py --sym-fc
+   % phono3py --fc-symmetry
    ```
 
-   where `--sym-fc` symmetrizes fc3 and fc2.
+   where `--fc-symmetry` symmetrizes fc3 and fc2.
 
 5) Calculate lattice thermal conductivity, e.g., by:
 
diff --git a/doc/random-displacements.md b/doc/random-displacements.md
index 4c4b1185..865d4c20 100644
--- a/doc/random-displacements.md
+++ b/doc/random-displacements.md
@@ -6,7 +6,7 @@ displacement-force dataset for computing force constants. This requires an
 external force constants calculator, e.g., symfc or ALM. Here, examples are
 presented with using symfc that can be installed via pip or conda easily.
 
-## Related setting tags
+## Related command options
 
 - {ref}`random_displacements_option` (`--rd`, `--random-seed`)
 - {ref}`fc_calculator_option` (`--fc-calc`)
@@ -82,27 +82,28 @@ calculated from `FORCES_FC3` (and optionally `FORCES_FC2`) and
 `phono3py_disp.yaml` by
 
 ```bash
-% phono3py --symfc -v
+% phono3py-load --symfc -v
 ```
 
-or with `phono3py_params.yaml`
+or
 
 ```bash
-% phono3py -c phono3py_params.yaml --symfc -v
+% phono3py-load phono3py_params.yaml --symfc -v
 ```
 
-Similarly, it is performed by also using `phono3py-load` command,
+Similarly, it is performed by also using `phono3py` command,
 
 ```bash
-% phono3py-load --symfc -v
+% phono3py --symfc -v
 ```
 
-or
+or with `phono3py_params.yaml`
 
 ```bash
-% phono3py-load phono3py_params.yaml --symfc -v
+% phono3py -c phono3py_params.yaml --symfc -v
 ```
 
+
 ## Cutoff pair-distance for fc3 calculation
 
 The number of supercells required for calculating fc3 depends on crystal
@@ -113,7 +114,7 @@ introducing cutoff distance for pairs of atoms. It is performed by
 `--fc-calc-opt` option as
 
 ```bash
-% phono3py --symfc -v --fc-calc-opt "cutoff=8"
+% phono3py-load --symfc -v --fc-calc-opt "cutoff=8"
 ```
 
 The shortcut of `--fc-calc-opt "cutoff=8"` is `--cutoff-pair 8`.
diff --git a/doc/turbomole.md b/doc/turbomole.md
index fe2af512..d007cf78 100644
--- a/doc/turbomole.md
+++ b/doc/turbomole.md
@@ -69,7 +69,7 @@ so the `-c control` parameter is not needed.
 4) Calculate 3rd and 2nd order force constants in files `fc3.hdf5` and `fc2.hdf5`:
 
    ```bash
-   % phono3py --turbomole --dim="2 2 2" --dim-fc2="3 3 3" --sym-fc
+   % phono3py --turbomole --dim="2 2 2" --dim-fc2="3 3 3" --fc-symmetry
    ```
 
    `--sym-fc` is used to symmetrize second- and third-order force constants.
diff --git a/doc/vasp.md b/doc/vasp.md
index be653e2d..e6cd650c 100644
--- a/doc/vasp.md
+++ b/doc/vasp.md
@@ -61,10 +61,10 @@
 4. Create `fc2.hdf` and `fc3.hdf`
 
    ```bash
-   % phono3py --sym-fc
+   % phono3py --fc-symmetry
    ```
 
-   `--sym-fc` symmetrizes fc3 and fc2. `fc2.hdf5` and `fc3.hdf5`
+   `--fc-symmetry` symmetrizes fc3 and fc2. `fc2.hdf5` and `fc3.hdf5`
    are created from `FORCES_FC3` (and
    optionally `FORCES_FC2`) and `phono3py_disp.yaml`. This step is
    not mandatory, but you can avoid calculating fc2 and fc3 at every
diff --git a/doc/wigner-solution.md b/doc/wigner-solution.md
index 234aeb22..ac6b1a93 100644
--- a/doc/wigner-solution.md
+++ b/doc/wigner-solution.md
@@ -38,7 +38,7 @@ As discussed in the references above, the term $\kappa_{\rm P}^{\alpha \beta}$ c
 To compute the Wigner conductivity with scattering in the RTA approximation, specify `--br` and `--wigner`. For `example/Wigner_La2Zr2O7`, the command is:
 
 ```bash
-phono3py --nac --cell POSCAR --fc2 --dim="2 2 2" --dim-fc2="4 4 4" --mesh="19 19 19" --tmin=300 --tmax=1000 --tstep=700 --sym-fc --isotope --br --wigner --read-gamma > tc_La2Zr2O7.out
+phono3py --nac --cell POSCAR --fc2 --dim="2 2 2" --dim-fc2="4 4 4" --mesh="19 19 19" --tmin=300 --tmax=1000 --tstep=700 --fc-symmetry --isotope --br --wigner --read-gamma > tc_La2Zr2O7.out
 ```
 The example above uses the `--read-gamma` option to read the phonon linewidths stored in the file  `kappa-m191919.hdf5`. The calculation of these linewidths is computationally expensive, more details are reported in the {ref}`paper on the Wigner formulation <citation_wigner_formulation>`, and in this example the linewidths are provided.
 To learn how to compute the linewidths, the reader is referred to the documentation of the `--write-gamma` option.
diff --git a/doc/workload-distribution.md b/doc/workload-distribution.md
index 4c7a8341..6d207f01 100644
--- a/doc/workload-distribution.md
+++ b/doc/workload-distribution.md
@@ -28,14 +28,14 @@ To avoid re-calculating fc3 and fc2, `fc3.hdf5` and `fc2.hdf5` are
 created on a single node:
 
 ```bash
-% phono3py --dim="2 2 2" --sym-fc -c POSCAR-unitcell
+% phono3py-load
 ```
 
 The indices of the irreducible grid-points neccesarry to specify
 `--ga` option are found by {ref}`--wgp option <wgp_option>`
 
 ```bash
-% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --wgp
+% phono3py-load --mesh 19 19 19 --br --wgp
 ```
 
 and they are stored in `ir_grid_points.yaml`.
@@ -51,7 +51,7 @@ computational demands into computer nodes, a set of the grid-point
 indices are chosen and executed as follows:
 
 ```bash
-% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --gp="0,1,2,3,4,5,6,7,8,9,20,21,22,23,24,25" --write-gamma
+% phono3py-load --mesh 19 19 19 --br --gp "0,1,2,3,4,5,6,7,8,9,20,21,22,23,24,25" --write-gamma
 ```
 
 Then many `kappa-m191919-gx.hdf5` files are generated. These file
@@ -63,7 +63,7 @@ indices, the RTA thermal conductivity is computed by another run in a
 short time from the stored data:
 
 ```bash
-% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --read-gamma
+% phono3py-load --mesh 19 19 19 --br --read-gamma
 ```
 
 ## A convenient script
@@ -96,7 +96,7 @@ with open("ir_grid_points.yaml") as f:
 Supposed that this script is saved as `divide_gps.py`,
 
 ```bash
-% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --wgp
+% phono3py-load --mesh 19 19 19 --wgp
 ...
 % python divide_gps.py 20
 0,30,52,82,120,402,434,468,524,844,1206
@@ -138,5 +138,5 @@ with `job.sh` (here for grid-engine):
 #$ -e err-phono3py-num.log
 #$ -o std-phono3py-num.log
 
-phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --gp="gps" --write-gamma
+phono3py-load --mesh 19 19 19 --br --gp "gps" --write-gamma
 ```
diff --git a/example/AlN-LDA/FORCES_FC3.lzma b/example/AlN-LDA/FORCES_FC3.lzma
deleted file mode 100644
index c9c35367..00000000
Binary files a/example/AlN-LDA/FORCES_FC3.lzma and /dev/null differ
diff --git a/example/AlN-LDA/FORCES_FC3.xz b/example/AlN-LDA/FORCES_FC3.xz
new file mode 100644
index 00000000..1f932734
Binary files /dev/null and b/example/AlN-LDA/FORCES_FC3.xz differ
diff --git a/example/AlN-LDA/README.md b/example/AlN-LDA/README.md
index 19dcdeb2..13e8cfd4 100644
--- a/example/AlN-LDA/README.md
+++ b/example/AlN-LDA/README.md
@@ -14,52 +14,56 @@ from all displaced supercell forces.
 
 Perfect and displaced supercells were created by
 
-```
+```bash
 % phono3py --dim 3 3 2 -c POSCAR-unitcell -d
 ```
 
-In the example directory, `FORCES_FC3` is compressed to `FORCES_FC3.lzma`. After
-unzipping `FORCES_FC3.lzma` (e.g., using `tar xvfz` or `tar xvfa`), to obtain
-`fc3.hdf5` and normal `fc2.hdf5`,
-
-```
-% phono3py --sym-fc
-```
-
-Using 13x13x9 sampling mesh, lattice thermal conductivity is calculated by
+In the example directory, `FORCES_FC3` is compressed to `FORCES_FC3.xz`. After
+unzipping `FORCES_FC3.xz` (e.g., using `xz -d`), to obtain `fc3.hdf5` and
+`fc2.hdf5` using symfc (the results without using symfc, i.e., finite difference
+method, are shown at the bottom of this README)
 
-```
-% phono3py --mesh 13 13 9 --fc3 --fc2 --br
+```bash
+% phono3py-load --symfc -v
 ```
 
-`kappa-m13139.hdf5` is written as the result. The lattice thermal conductivity
-is calculated as k_xx=228.2 and k_zz=224.1 W/m-K at 300 K.
+Lattice thermal conductivity is calculated by
 
-With `--nac` option, non-analytical term correction is applied reading the Born
-effective charges and dielectric constant from `BORN` file:
-
-```
-% phono3py --mesh 13 13 9 --fc3 --fc2 --br --nac
+```bash
+% phono3py-load --mesh 40 --br --ts 300
 ```
 
-This changes thermal conductivity at 300 K to k_xx=235.7 and k_zz=219.1. The
-shape of phonon band structure is important to fullfil energy and momentum
-conservations.
+`kappa-m15158.hdf5` is written as the result. Parameters for non-analytical term
+correction (NAC) is automatically read from those stored in `phono3py_disp.yaml` or
+`BORN` file. The lattice thermal conductivity is calculated as k_xx=242.8 and
+k_zz=226.5 W/m-K at 300 K. Without NAC, k_xx=233.6 and k_zz=222.2.
 
-Use of larger supercell of fc2 may change the shape of phonon band structure. To
-see it, first regenerate `phono3py_disp.yaml` with `--dim-fc2` option,
+Use of larger supercell for fc2 may change the shape of phonon band structure.
+To see it, first regenerate `phono3py_disp.yaml` with `--dim-fc2` option,
 
-```
+```bash
 % phono3py --dim 3 3 2 --dim-fc2 5 5 3 -c POSCAR-unitcell -d
 ```
 
 Then re-create force constants and calculate thermal conductivity,
 
+```bash
+% phono3py-load --symfc -v
+% phono3py-load --br --mesh=40 --ts 300
 ```
-% phono3py --sym-fc
-% phono3py --mesh="13 13 9" --fc3 --fc2 --br --nac
+
+If `phono3py_disp.yaml` is renamed to `phono3py_disp_dimfc2.yaml`, it can be
+specified at the first argument of `phono3py-load` command:
+
+```bash
+% phono3py-load phono3py_disp_dimfc2.yaml --symfc -v
+% phono3py-load phono3py_disp_dimfc2.yaml --br --mesh=40 --ts 300
 ```
 
-k_xx=236.0 and k_zz=222.2 are obtained. In the case of this example, we can see
+k_xx=240.2 and k_zz=230.1 are obtained. In the case of this example, we can see
 that the larger fc2 supercell contributes little, which means that the 3x3x2
 supercell was good enough to obtain a good shape of phonon band structure.
+
+Using the finite difference method implemented in phono3py, lattice thermal
+conductivities are obtained as k_xx=251.2 and k_zz=233,4 without using the large
+fc2 supercell and k_xx=249.4 k_zz=236.9 using the large fc2 supercell.
diff --git a/example/NaCl-alm/README.md b/example/NaCl-alm/README.md
index bc0bf93b..1a75f98e 100644
--- a/example/NaCl-alm/README.md
+++ b/example/NaCl-alm/README.md
@@ -29,4 +29,4 @@ lattice thermal conductivity at 300 K is calculated by
 % phono3py-load phono3py_params_NaCl222.yaml.xz --mesh 50 --ts 300 --br
 ```
 
-The result is ~7.2 W/m-K at 300 K.
+The result is ~7.3 W/m-K at 300 K.
diff --git a/example/Si-CRYSTAL/README b/example/Si-CRYSTAL/README
index d8f15907..e31a04d5 100644
--- a/example/Si-CRYSTAL/README
+++ b/example/Si-CRYSTAL/README
@@ -23,7 +23,7 @@ for the CRYSTAL interface, so the -c crystal.o parameter is not needed
 
 4) Create force constant files fc2.hdf5 and fc3.hdf5:
 
-   phono3py --crystal --dim="2 2 2" --dim-fc2="4 4 4" --sym-fc
+   phono3py --crystal --dim="2 2 2" --dim-fc2="4 4 4" --fc-symmetry
 
 5) Thermal conductivity calculation:
 
diff --git a/example/Si-LDA/README.md b/example/Si-LDA/README.md
index 751ddffd..8891834b 100644
--- a/example/Si-LDA/README.md
+++ b/example/Si-LDA/README.md
@@ -13,18 +13,18 @@ conventional unit cell to the primitive cell.
 To create `fc3.hdf5` and `fc2.hdf5`,
 
 ```
-% phono3py --sym-fc
+% phono3py-load
 ```
 
 Using 11x11x11 sampling mesh, lattice thermal conductivity is calculated by
 
 ```
-% phono3py --mesh 11 11 11 --fc3 --fc2 --br
+% phono3py-load --mesh 11 11 11 --br --ts 300
 ```
 
 `kappa-m111111.hdf5` is written as the result. The lattice thermal conductivity
-is calculated as 112.5 W/m-K at 300 K. This becomes, with 19x19x19 sampling
-mesh, 127.0 W/m-K.
+is calculated as 112.4 W/m-K at 300 K. This becomes, with 19x19x19 sampling
+mesh, 128.2 W/m-K.
 
 Accumulated lattice thermal conductivity is calculated using `phono3py-kaccum`
 script.
diff --git a/example/Si-PBE/README.md b/example/Si-PBE/README.md
index 89cd4195..a9b62b79 100644
--- a/example/Si-PBE/README.md
+++ b/example/Si-PBE/README.md
@@ -13,18 +13,18 @@ conventional unit cell to the primitive cell.
 To create `fc3.hdf5` and `fc2.hdf5`,
 
 ```
-% phono3py --sym-fc
+% phono3py-load
 ```
 
 Using 11x11x11 sampling mesh, lattice thermal conductivity is calculated by
 
 ```
-% phono3py --mesh 11 11 11 --fc3 --fc2 --br
+% phono3py-load --mesh 11 11 11 --fc3 --fc2 --br
 ```
 
 `kappa-m111111.hdf5` is written as the result. The lattice thermal conductivity
-is calculated as 119.3 W/m-K at 300 K. This becomes, with 19x19x19 sampling
-mesh, 132.4 W/m-K.
+is calculated as 119.5 W/m-K at 300 K. This becomes, with 19x19x19 sampling
+mesh, 130.1 W/m-K.
 
 Accumulated lattice thermal conductivity is calculated using `phono3py-kaccum`
 script.
diff --git a/example/Si-PBEsol/README.md b/example/Si-PBEsol/README.md
index aed2fe8c..ee16bfba 100644
--- a/example/Si-PBEsol/README.md
+++ b/example/Si-PBEsol/README.md
@@ -13,18 +13,18 @@ the conventional unit cell to the primitive cell.
 To create `fc3.hdf5` and `fc2.hdf5`,
 
 ```
-% phono3py --sym-fc
+% phono3py-load
 ```
 
 Using 11x11x11 sampling mesh, lattice thermal conductivity is calculated by
 
 ```
-% phono3py --mesh 11 11 11 --fc3 --fc2 --br
+% phono3py-load --mesh 11 11 11 --br
 ```
 
 `kappa-m111111.hdf5` is written as the result. The lattice thermal conductivity
-is calculated as 109.0 W/m-K at 300 K. This becomes, with 19x19x19 sampling
-mesh, 123.2 W/m-K.
+is calculated as 109.1 W/m-K at 300 K. This becomes, with 19x19x19 sampling
+mesh, 124.4 W/m-K.
 
 Accumulated lattice thermal conductivity is calculated using `phono3py-kaccum`
 script.
@@ -80,11 +80,11 @@ In this example directory, the dataset is renamed to
 at `vasprun_xml_fc2/disp-fc2-00001/vasprun.xml`. `FORCES_FC2` is generated by
 
 ```
-% phono3py phono3py_disp_dimfc2.yaml --cf2 vasprun_xml_fc2/disp-fc2-00001/vasprun.xml
+% phono3py -c phono3py_disp_dimfc2.yaml --cf2 vasprun_xml_fc2/disp-fc2-00001/vasprun.xml
 ```
 
 A lattice thermal conductivity calculation is performed by, e.g.,
 ```
 % phono3py-load phono3py_disp_dimfc2.yaml --mesh 11 11 11 --br --ts 300
 ```
-The result is 108.0 W/m-K, and with the 19x19x19 mesh, 125.4 W/m-K.
+The result is 107.9 W/m-K, and with the 19x19x19 mesh, 125.4 W/m-K.
diff --git a/example/Si-QE/README.md b/example/Si-QE/README.md
index 8f17a517..1b8d5b88 100644
--- a/example/Si-QE/README.md
+++ b/example/Si-QE/README.md
@@ -13,18 +13,18 @@ the conventional unit cell to the primitive cell.
 To create `fc3.hdf5` and `fc2.hdf5`,
 
 ```
-% phono3py --sym-fc
+% phono3py-load
 ```
 
 Using 11x11x11 sampling mesh, lattice thermal conductivity is calculated by
 
 ```
-% phono3py --mesh 11 11 11 --fc3 --fc2 --br
+% phono3py-load --mesh 11 11 11 --br --ts 300
 ```
 
 `kappa-m111111.hdf5` is written as the result. The lattice thermal conductivity
-is calculated as 118.9 W/m-K at 300 K. This becomes, with 19x19x19 sampling
-mesh, 129.9 W/m-K.
+is calculated as 121.3 W/m-K at 300 K. This becomes, with 19x19x19 sampling
+mesh, 135.4 W/m-K.
 
 The .out files for supercells are found in `supercell_out.tar.lzma`. If phono3py
 is properly installed, the following command should work.
diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py
index 294ac2dd..00d54efe 100644
--- a/phono3py/api_phono3py.py
+++ b/phono3py/api_phono3py.py
@@ -1487,6 +1487,9 @@ def produce_fc3(
             fc2 = None
         if fc_calculator_options is not None and "|" in fc_calculator_options:
             fc2 = None
+        # fc2 should not be set if phonon_supercell_matrix is available.
+        if self._phonon_supercell_matrix is not None:
+            fc2 = None
 
         # Normally self._fc2 is overwritten in produce_fc2
         if self._fc2 is None:
@@ -1527,6 +1530,9 @@ def produce_fc2(
         else:
             disp_dataset = self._phonon_dataset
 
+        if not forces_in_dataset(disp_dataset):
+            raise RuntimeError("Forces are not set in the dataset.")
+
         if is_compact_fc:
             p2s_map = self._phonon_primitive.p2s_map
         else:
@@ -1548,10 +1554,7 @@ def produce_fc2(
             )
         else:
             if "displacements" in disp_dataset:
-                msg = (
-                    "fc_calculator has to be set to produce force "
-                    "constans from this dataset for fc2."
-                )
+                msg = "fc_calculator to solve fc2 has to be set."
                 raise RuntimeError(msg)
             self._fc2 = get_phonopy_fc2(
                 self._phonon_supercell,
diff --git a/phono3py/cui/create_force_constants.py b/phono3py/cui/create_force_constants.py
index 740455f0..8ac34d9f 100644
--- a/phono3py/cui/create_force_constants.py
+++ b/phono3py/cui/create_force_constants.py
@@ -37,11 +37,10 @@
 from __future__ import annotations
 
 import copy
-import os
 import pathlib
 import sys
 from dataclasses import asdict
-from typing import Optional, Union
+from typing import Literal, Optional, Union
 
 import numpy as np
 from phonopy.cui.phonopy_script import file_exists, print_error
@@ -81,6 +80,7 @@ def create_phono3py_force_constants(
     settings,
     ph3py_yaml: Optional[Phono3pyYaml] = None,
     phono3py_yaml_filename: Optional[str] = None,
+    calculator: Optional[str] = None,
     input_filename: Optional[str] = None,
     output_filename: Optional[str] = None,
     log_level=1,
@@ -118,15 +118,12 @@ def create_phono3py_force_constants(
         if settings.read_fc3:
             _read_phono3py_fc3(phono3py, symmetrize_fc3r, input_filename, log_level)
         else:  # fc3 from FORCES_FC3 or ph3py_yaml
-            _create_phono3py_fc3(
+            _read_dataset_fc3(
                 phono3py,
                 ph3py_yaml,
                 phono3py_yaml_filename,
-                symmetrize_fc3r,
-                settings.is_compact_fc,
                 settings.cutoff_pair_distance,
-                fc_calculator,
-                fc_calculator_options,
+                calculator,
                 settings.use_pypolymlp,
                 settings.mlp_params,
                 settings.displacement_distance,
@@ -134,6 +131,14 @@ def create_phono3py_force_constants(
                 settings.random_seed,
                 log_level,
             )
+            phono3py.produce_fc3(
+                symmetrize_fc3r=symmetrize_fc3r,
+                is_compact_fc=settings.is_compact_fc,
+                fc_calculator=extract_fc2_fc3_calculators(fc_calculator, 3),
+                fc_calculator_options=extract_fc2_fc3_calculators(
+                    fc_calculator_options, 3
+                ),
+            )
 
         cutoff_distance = settings.cutoff_fc3_distance
         if cutoff_distance is not None and cutoff_distance > 0:
@@ -168,15 +173,20 @@ def create_phono3py_force_constants(
     if settings.read_fc2:
         _read_phono3py_fc2(phono3py, symmetrize_fc2, input_filename, log_level)
     else:
-        _create_phono3py_fc2(
-            phono3py,
-            ph3py_yaml,
-            symmetrize_fc2,
-            settings.is_compact_fc,
-            fc_calculator,
-            fc_calculator_options,
-            log_level,
+        if phono3py.dataset is None or phono3py.phonon_supercell_matrix is not None:
+            _read_dataset_fc2(
+                phono3py,
+                ph3py_yaml,
+                calculator,
+                log_level,
+            )
+        phono3py.produce_fc2(
+            symmetrize_fc2=symmetrize_fc2,
+            is_compact_fc=settings.is_compact_fc,
+            fc_calculator=extract_fc2_fc3_calculators(fc_calculator, 2),
+            fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 2),
         )
+
         if output_filename is None:
             filename = "fc2.hdf5"
         else:
@@ -201,7 +211,8 @@ def parse_forces(
     cutoff_pair_distance=None,
     force_filename: str = "FORCES_FC3",
     phono3py_yaml_filename: Optional[str] = None,
-    fc_type=None,
+    fc_type: Literal["fc3", "phonon_fc2"] = "fc3",
+    calculator: Optional[str] = None,
     log_level=0,
 ):
     """Read displacements and forces.
@@ -215,53 +226,39 @@ def parse_forces(
     """
     filename_read_from: Optional[str] = None
     dataset = None
-    calculator = phono3py.calculator
+
+    if phono3py.phonon_supercell is None or fc_type == "fc3":
+        natom = len(phono3py.supercell)
+    else:
+        natom = len(phono3py.phonon_supercell)
 
     # Get dataset from ph3py_yaml. dataset can be None.
     # physical_units can be overwritten if calculator is found in ph3py_yaml.
     if ph3py_yaml:
         dataset = _extract_dataset_from_ph3py_yaml(ph3py_yaml, fc_type)
-    if dataset and ph3py_yaml.calculator:
-        calculator = ph3py_yaml.calculator
+        if dataset:
+            filename_read_from = phono3py_yaml_filename
 
     physical_units = get_default_physical_units(calculator)
 
-    if phono3py.phonon_supercell is None or fc_type == "fc3":
-        natom = len(phono3py.supercell)
-    else:
-        natom = len(phono3py.phonon_supercell)
-
-    if dataset:
-        filename_read_from = phono3py_yaml_filename
-
-        # Units of displacements and forces are converted. If forces don't
-        # exist, the convesion will not be performed for forces.
-        if calculator is not None:
-            _convert_unit_in_dataset(
-                dataset,
-                distance_to_A=physical_units["distance_to_A"],
-                force_to_eVperA=physical_units["force_to_eVperA"],
+    # Forces are not yet found in dataset. Then try to read from FORCES_FC3 or
+    # FORCES_FC2.
+    if force_filename is not None:
+        if dataset is None or (dataset is not None and not forces_in_dataset(dataset)):
+            dataset = _read_FORCES_FC3_or_FC2(
+                natom, dataset, fc_type, filename=force_filename, log_level=log_level
             )
-
-    # Try to read FORCES_FC* if type-2 and return dataset.
-    # None is returned unless type-2.
-    # can emit FileNotFoundError.
-    if dataset is None or (dataset is not None and not forces_in_dataset(dataset)):
-        _dataset = read_type2_dataset(
-            natom, filename=force_filename, log_level=log_level
+            if dataset:
+                filename_read_from = force_filename
+
+    # Units of displacements and forces are converted. If forces don't
+    # exist, the convesion will not be performed for forces.
+    if calculator is not None:
+        _convert_unit_in_dataset(
+            dataset,
+            distance_to_A=physical_units["distance_to_A"],
+            force_to_eVperA=physical_units["force_to_eVperA"],
         )
-        # Do not overwrite dataset when _dataset is None.
-        if _dataset:
-            filename_read_from = force_filename
-            dataset = _dataset
-
-            # Units of displacements and forces are converted.
-            if calculator is not None:
-                _convert_unit_in_dataset(
-                    dataset,
-                    distance_to_A=physical_units["distance_to_A"],
-                    force_to_eVperA=physical_units["force_to_eVperA"],
-                )
 
     assert dataset is not None
 
@@ -286,29 +283,6 @@ def parse_forces(
             if log_level:
                 print("Cutoff-pair-distance: %f" % cutoff_pair_distance)
 
-    # Type-1 FORCES_FC*.
-    # dataset comes either from disp_fc*.yaml or phono3py*.yaml.
-    if not forces_in_dataset(dataset):
-        if force_filename is not None:
-            if fc_type == "fc3":
-                parse_FORCES_FC3(dataset, filename=force_filename)
-            else:
-                parse_FORCES_FC2(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.
-        if calculator is not None:
-            _convert_unit_in_dataset(
-                dataset,
-                force_to_eVperA=physical_units["force_to_eVperA"],
-            )
-
     return dataset
 
 
@@ -437,32 +411,49 @@ def _read_phono3py_fc2(phono3py, symmetrize_fc2, input_filename, log_level):
     phono3py.fc2 = phonon_fc2
 
 
-def read_type2_dataset(natom, filename="FORCES_FC3", log_level=0) -> Optional[dict]:
-    """Read type-2 FORCES_FC3."""
+def _read_FORCES_FC3_or_FC2(
+    natom: int,
+    dataset: Optional[dict],
+    fc_type: str,
+    filename: str = "FORCES_FC3",
+    log_level: int = 0,
+) -> Optional[dict]:
+    """Read FORCES_FC3 or FORCES_FC2.
+
+    Read the first line of forces file to determine the type of the file.
+
+    """
     if filename is None or not pathlib.Path(filename).exists():
         return None
 
     with open(filename, "r") as f:
         len_first_line = get_length_of_first_line(f)
-        if len_first_line == 6:
-            dataset = get_dataset_type2(f, natom)
+        if len_first_line == 6:  # Type-2
+            _dataset = get_dataset_type2(f, natom)
             if log_level:
-                n_disp = len(dataset["displacements"])
+                n_disp = len(_dataset["displacements"])
                 print(f'{n_disp} snapshots were found in "{filename}".')
-        else:
-            dataset = None
+            return _dataset
+
+    # Type-1
+    if fc_type == "fc3":
+        parse_FORCES_FC3(dataset, filename)
+    else:
+        parse_FORCES_FC2(dataset, filename)
+    if log_level:
+        print(
+            f'Sets of supercell forces were read from "{filename}".',
+            flush=True,
+        )
     return dataset
 
 
-def _create_phono3py_fc3(
+def _read_dataset_fc3(
     phono3py: Phono3py,
     ph3py_yaml: Optional[Phono3pyYaml],
     phono3py_yaml_filename: Optional[str],
-    symmetrize_fc3r: bool,
-    is_compact_fc: bool,
     cutoff_pair_distance: Optional[float],
-    fc_calculator: Optional[str],
-    fc_calculator_options: Optional[str],
+    calculator: Optional[str],
     use_pypolymlp: bool,
     mlp_params: Union[str, dict, PypolymlpParams],
     displacement_distance: Optional[float],
@@ -484,16 +475,15 @@ def _create_phono3py_fc3(
     when the former value is smaller than the later.
 
     """
-    _ph3py_yaml = _get_default_ph3py_yaml(ph3py_yaml)
-
     try:
         dataset = parse_forces(
             phono3py,
-            ph3py_yaml=_ph3py_yaml,
+            ph3py_yaml=ph3py_yaml,
             cutoff_pair_distance=cutoff_pair_distance,
             force_filename="FORCES_FC3",
             phono3py_yaml_filename=phono3py_yaml_filename,
             fc_type="fc3",
+            calculator=calculator,
             log_level=log_level,
         )
     except RuntimeError as e:
@@ -518,12 +508,6 @@ def _create_phono3py_fc3(
         )
     else:
         phono3py.dataset = dataset
-    phono3py.produce_fc3(
-        symmetrize_fc3r=symmetrize_fc3r,
-        is_compact_fc=is_compact_fc,
-        fc_calculator=extract_fc2_fc3_calculators(fc_calculator, 3),
-        fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 3),
-    )
 
 
 def run_pypolymlp_to_compute_forces(
@@ -653,13 +637,10 @@ def run_pypolymlp_to_compute_phonon_forces(
     ph3py.evaluate_phonon_mlp()
 
 
-def _create_phono3py_fc2(
+def _read_dataset_fc2(
     phono3py: Phono3py,
     ph3py_yaml: Optional[Phono3pyYaml],
-    symmetrize_fc2,
-    is_compact_fc,
-    fc_calculator,
-    fc_calculator_options,
+    calculator,
     log_level,
 ):
     """Read forces and produce fc2.
@@ -667,19 +648,22 @@ def _create_phono3py_fc2(
     force_filename is either "FORCES_FC2" or "FORCES_FC3".
 
     """
-    if phono3py.phonon_supercell_matrix is None:
+    # _ph3py_yaml = _get_default_ph3py_yaml(ph3py_yaml)
+
+    if phono3py.phonon_supercell_matrix is not None:
+        force_filename = "FORCES_FC2"
+    elif phono3py.dataset is None:
         force_filename = "FORCES_FC3"
     else:
-        force_filename = "FORCES_FC2"
-
-    _ph3py_yaml = _get_default_ph3py_yaml(ph3py_yaml)
+        raise RuntimeError("Force filename is not determined.")
 
     try:
         dataset = parse_forces(
             phono3py,
-            ph3py_yaml=_ph3py_yaml,
+            ph3py_yaml=ph3py_yaml,
             force_filename=force_filename,
             fc_type="phonon_fc2",
+            calculator=calculator,
             log_level=log_level,
         )
     except RuntimeError as e:
@@ -690,18 +674,15 @@ def _create_phono3py_fc2(
     except FileNotFoundError as e:
         file_exists(e.filename, log_level)
 
-    phono3py.phonon_dataset = dataset
-    phono3py.produce_fc2(
-        symmetrize_fc2=symmetrize_fc2,
-        is_compact_fc=is_compact_fc,
-        fc_calculator=extract_fc2_fc3_calculators(fc_calculator, 2),
-        fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 2),
-    )
+    if phono3py.phonon_supercell_matrix is not None:
+        phono3py.phonon_dataset = dataset
+    elif phono3py.dataset is None:
+        phono3py.dataset = dataset
 
 
 def _get_default_ph3py_yaml(ph3py_yaml: Optional[Phono3pyYaml]):
     _ph3py_yaml = ph3py_yaml
-    if _ph3py_yaml is None and os.path.isfile("phono3py_disp.yaml"):
+    if _ph3py_yaml is None and pathlib.Path("phono3py_disp.yaml").exists():
         _ph3py_yaml = Phono3pyYaml()
         _ph3py_yaml.read("phono3py_disp.yaml")
     return _ph3py_yaml
diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py
index d100afdb..f87d6e2c 100644
--- a/phono3py/cui/load.py
+++ b/phono3py/cui/load.py
@@ -46,6 +46,7 @@
 from phonopy.interface.calculator import get_default_physical_units
 from phonopy.structure.atoms import PhonopyAtoms
 from phonopy.structure.cells import determinant
+from phonopy.units import VaspToTHz
 
 from phono3py import Phono3py
 from phono3py.cui.create_force_constants import (
@@ -261,6 +262,7 @@ def load(
         or unitcell is not None
         or unitcell_filename is not None
     ):
+        _calculator = calculator
         cell, smat, pmat = load_helper.get_cell_settings(
             supercell_matrix=supercell_matrix,
             primitive_matrix=primitive_matrix,
@@ -268,7 +270,7 @@ def load(
             supercell=supercell,
             unitcell_filename=unitcell_filename,
             supercell_filename=supercell_filename,
-            calculator=calculator,
+            calculator=_calculator,
             symprec=symprec,
         )
         if phonon_supercell_matrix is not None:
@@ -286,7 +288,8 @@ def load(
     elif phono3py_yaml is not None:
         ph3py_yaml = Phono3pyYaml()
         ph3py_yaml.read(phono3py_yaml)
-        cell = ph3py_yaml.unitcell
+        cell = ph3py_yaml.unitcell.copy()
+        _calculator = ph3py_yaml.calculator
         smat = ph3py_yaml.supercell_matrix
         ph_smat = ph3py_yaml.phonon_supercell_matrix
         if smat is None:
@@ -303,10 +306,13 @@ def load(
         else:
             _nac_params = None
 
-    # units keywords: factor, nac_factor, distance_to_A
-    physical_units = get_default_physical_units(calculator)
+    # Convert distance unit of unit cell to Angstrom
+    physical_units = get_default_physical_units(_calculator)
+    factor_to_A = physical_units["distance_to_A"]
+    cell.cell = cell.cell * factor_to_A
+
     if factor is None:
-        _factor = physical_units["factor"]
+        _factor = VaspToTHz
     else:
         _factor = factor
     ph3py = Phono3py(
@@ -320,7 +326,6 @@ def load(
         is_mesh_symmetry=is_mesh_symmetry,
         use_grg=use_grg,
         make_r0_average=make_r0_average,
-        calculator=calculator,
         log_level=log_level,
     )
 
@@ -343,6 +348,7 @@ def load(
         forces_fc3_filename=forces_fc3_filename,
         forces_fc2_filename=forces_fc2_filename,
         phono3py_yaml_filename=phono3py_yaml,
+        calculator=_calculator,
         use_pypolymlp=use_pypolymlp,
         log_level=log_level,
     )
@@ -379,6 +385,7 @@ def set_dataset_and_force_constants(
     forces_fc2_filename: Optional[Union[os.PathLike, Sequence]] = None,
     phono3py_yaml_filename: Optional[os.PathLike] = None,
     cutoff_pair_distance: Optional[float] = None,
+    calculator: Optional[str] = None,
     use_pypolymlp: bool = False,
     log_level: int = 0,
 ) -> dict:
@@ -404,6 +411,7 @@ def set_dataset_and_force_constants(
         forces_fc3_filename=forces_fc3_filename,
         phono3py_yaml_filename=phono3py_yaml_filename,
         cutoff_pair_distance=cutoff_pair_distance,
+        calculator=calculator,
         log_level=log_level,
     )
     if not read_fc["fc3"]:
@@ -416,6 +424,7 @@ def set_dataset_and_force_constants(
         ph3py_yaml=ph3py_yaml,
         fc2_filename=fc2_filename,
         forces_fc2_filename=forces_fc2_filename,
+        calculator=calculator,
         log_level=log_level,
     )
     if not read_fc["fc2"]:
@@ -475,17 +484,23 @@ def compute_force_constants_from_datasets(
             if log_level and symmetrize_fc and fc_calculator is None:
                 print("fc3 was symmetrized.")
 
-    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,
-            fc_calculator=fc2_calculator,
-            fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 2),
-        )
-        if log_level and symmetrize_fc and fc_calculator is None:
-            print("fc2 was symmetrized.")
+    if not read_fc["fc2"]:
+        if (
+            ph3py.phonon_supercell_matrix is None and forces_in_dataset(ph3py.dataset)
+        ) or (
+            ph3py.phonon_supercell_matrix is not None
+            and forces_in_dataset(ph3py.phonon_dataset)
+        ):
+            ph3py.produce_fc2(
+                symmetrize_fc2=symmetrize_fc,
+                is_compact_fc=is_compact_fc,
+                fc_calculator=fc2_calculator,
+                fc_calculator_options=extract_fc2_fc3_calculators(
+                    fc_calculator_options, 2
+                ),
+            )
+            if log_level and symmetrize_fc and fc_calculator is None:
+                print("fc2 was symmetrized.")
 
 
 def _get_dataset_or_fc3(
@@ -495,6 +510,7 @@ def _get_dataset_or_fc3(
     forces_fc3_filename: Optional[Union[os.PathLike, Sequence]] = None,
     phono3py_yaml_filename: Optional[os.PathLike] = None,
     cutoff_pair_distance: Optional[float] = None,
+    calculator: Optional[str] = None,
     log_level: int = 0,
 ) -> tuple[bool, dict]:
     p2s_map = ph3py.primitive.p2s_map
@@ -523,6 +539,7 @@ def _get_dataset_or_fc3(
             None,
             phono3py_yaml_filename,
             cutoff_pair_distance,
+            calculator,
             log_level,
         )
     elif forces_fc3_filename is not None or pathlib.Path("FORCES_FC3").exists():
@@ -536,6 +553,7 @@ def _get_dataset_or_fc3(
             force_filename,
             phono3py_yaml_filename,
             cutoff_pair_distance,
+            calculator,
             log_level,
         )
     elif ph3py_yaml is not None and ph3py_yaml.dataset is not None:
@@ -547,6 +565,7 @@ def _get_dataset_or_fc3(
             None,
             phono3py_yaml_filename,
             cutoff_pair_distance,
+            calculator,
             log_level,
         )
 
@@ -558,6 +577,7 @@ def _get_dataset_phonon_dataset_or_fc2(
     ph3py_yaml: Optional[Phono3pyYaml] = None,
     fc2_filename: Optional[os.PathLike] = None,
     forces_fc2_filename: Optional[Union[os.PathLike, Sequence]] = None,
+    calculator: Optional[str] = None,
     log_level: int = 0,
 ) -> tuple[bool, dict, dict]:
     phonon_p2s_map = ph3py.phonon_primitive.p2s_map
@@ -584,6 +604,7 @@ def _get_dataset_phonon_dataset_or_fc2(
             ph3py_yaml,
             None,
             "phonon_fc2",
+            calculator,
             log_level,
         )
     elif (
@@ -598,6 +619,7 @@ def _get_dataset_phonon_dataset_or_fc2(
             ph3py_yaml,
             force_filename,
             "phonon_fc2",
+            calculator,
             log_level,
         )
     elif ph3py_yaml is not None and ph3py_yaml.phonon_dataset is not None:
@@ -608,6 +630,7 @@ def _get_dataset_phonon_dataset_or_fc2(
             ph3py_yaml,
             None,
             "phonon_fc2",
+            calculator,
             log_level,
         )
 
@@ -620,6 +643,7 @@ def _get_dataset_for_fc3(
     force_filename,
     phono3py_yaml_filename,
     cutoff_pair_distance,
+    calculator,
     log_level,
 ) -> dict:
     dataset = parse_forces(
@@ -629,6 +653,7 @@ def _get_dataset_for_fc3(
         force_filename=force_filename,
         phono3py_yaml_filename=phono3py_yaml_filename,
         fc_type="fc3",
+        calculator=calculator,
         log_level=log_level,
     )
     return dataset
@@ -639,6 +664,7 @@ def _get_dataset_for_fc2(
     ph3py_yaml: Optional[Phono3pyYaml],
     force_filename,
     fc_type,
+    calculator,
     log_level,
 ):
     dataset = parse_forces(
@@ -646,6 +672,7 @@ def _get_dataset_for_fc2(
         ph3py_yaml=ph3py_yaml,
         force_filename=force_filename,
         fc_type=fc_type,
+        calculator=calculator,
         log_level=log_level,
     )
     return dataset
diff --git a/phono3py/cui/phono3py_argparse.py b/phono3py/cui/phono3py_argparse.py
index 0c31dde6..32c392f4 100644
--- a/phono3py/cui/phono3py_argparse.py
+++ b/phono3py/cui/phono3py_argparse.py
@@ -59,13 +59,14 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
         default=False,
         help=("Use ALM for generating 2nd and 3rd force constants " "in one fitting"),
     )
-    parser.add_argument(
-        "--amplitude",
-        dest="displacement_distance",
-        type=float,
-        default=None,
-        help="Distance of displacements",
-    )
+    if not load_phono3py_yaml:
+        parser.add_argument(
+            "--amplitude",
+            dest="displacement_distance",
+            type=float,
+            default=None,
+            help="Distance of displacements",
+        )
     parser.add_argument(
         "--ave-pp",
         dest="use_ave_pp",
diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py
index 0b0c6c2c..b0bc4821 100644
--- a/phono3py/cui/phono3py_script.py
+++ b/phono3py/cui/phono3py_script.py
@@ -38,6 +38,7 @@
 
 import argparse
 import datetime
+import pathlib
 import sys
 from typing import Optional
 
@@ -92,6 +93,7 @@
 )
 from phono3py.interface.phono3py_yaml import Phono3pyYaml
 from phono3py.phonon.grid import get_grid_point_from_address, get_ir_grid_points
+from phono3py.phonon3.dataset import forces_in_dataset
 from phono3py.phonon3.fc3 import show_drift_fc3
 from phono3py.phonon3.gruneisen import run_gruneisen_parameters
 from phono3py.version import __version__
@@ -485,7 +487,6 @@ def init_phono3py(
         use_grg=settings.use_grg,
         make_r0_average=settings.is_fc3_r0_average,
         symprec=symprec,
-        calculator=interface_mode,
         log_level=log_level,
     )
     phono3py.masses = settings.masses
@@ -541,15 +542,16 @@ def create_supercells_with_displacements(
             log_level=log_level,
         )
 
-        store_nac_params(
-            phono3py,
-            settings,
-            cell_info["phonopy_yaml"],
-            unitcell_filename,
-            log_level,
-            nac_factor=Hartree * Bohr,
-            load_phonopy_yaml=load_phono3py_yaml,
-        )
+        if pathlib.Path("BORN").exists():
+            store_nac_params(
+                phono3py,
+                settings,
+                cell_info["phonopy_yaml"],
+                unitcell_filename,
+                log_level,
+                nac_factor=Hartree * Bohr,
+                load_phonopy_yaml=load_phono3py_yaml,
+            )
 
         if log_level:
             if phono3py.supercell.magnetic_moments is None:
@@ -574,6 +576,7 @@ def store_force_constants(
     settings,
     ph3py_yaml: Phono3pyYaml,
     phono3py_yaml_filename,
+    calculator,
     input_filename,
     output_filename,
     load_phono3py_yaml,
@@ -590,6 +593,7 @@ def store_force_constants(
             phono3py_yaml_filename=phono3py_yaml_filename,
             cutoff_pair_distance=settings.cutoff_pair_distance,
             use_pypolymlp=settings.use_pypolymlp,
+            calculator=calculator,
             log_level=log_level,
         )
         (fc_calculator, fc_calculator_options) = get_fc_calculator_params(
@@ -616,10 +620,18 @@ def store_force_constants(
         if log_level:
             if phono3py.fc3 is None:
                 print("fc3 could not be obtained.")
+                if not forces_in_dataset(phono3py.dataset):
+                    print("Forces were not found.")
             else:
                 show_drift_fc3(phono3py.fc3, primitive=phono3py.primitive)
             if phono3py.fc2 is None:
                 print("fc2 could not be obtained.")
+                if phono3py.phonon_supercell_matrix is None:
+                    if not forces_in_dataset(phono3py.dataset):
+                        print("Forces were not found.")
+                else:
+                    if not forces_in_dataset(phono3py.phonon_dataset):
+                        print("Forces for dim-fc2 were not found.")
             else:
                 show_drift_force_constants(
                     phono3py.fc2, primitive=phono3py.phonon_primitive, name="fc2"
@@ -661,6 +673,7 @@ def store_force_constants(
                 settings,
                 ph3py_yaml=ph3py_yaml,
                 phono3py_yaml_filename=phono3py_yaml_filename,
+                calculator=calculator,
                 input_filename=input_filename,
                 output_filename=output_filename,
                 log_level=log_level,
@@ -1129,6 +1142,7 @@ def main(**argparse_control):
         settings,
         cell_info["phonopy_yaml"],
         unitcell_filename,
+        interface_mode,
         input_filename,
         output_filename,
         load_phono3py_yaml,
diff --git a/phono3py/interface/phono3py_yaml.py b/phono3py/interface/phono3py_yaml.py
index 36fdff8e..cc6cc13b 100644
--- a/phono3py/interface/phono3py_yaml.py
+++ b/phono3py/interface/phono3py_yaml.py
@@ -133,10 +133,34 @@ def _parse_dataset(self):
 
         This method override PhonopyYaml._parse_dataset.
 
+        Phonon_displacements in type1 in old format is represented in yaml by
+
+        displacements:
+
+        But displacements (for fc3) in type2 in old format is also represented
+        by
+
+        displacements:
+
+        Therefore, these have to be distinguished by phonon_supercell_matrix.
+
+        In new format of type1
+
+        displacements: -> phonon_displacements:
+
+        and displacements: of type2 is put under phonon_dataset: block.
+
         """
-        self._data.phonon_dataset = self._get_dataset(
-            self._data.phonon_supercell, key_prefix="phonon_"
-        )
+        if (
+            self._data.phonon_supercell_matrix is not None
+            and "phonon_displacements" not in self._yaml
+            and "displacements" in self._yaml
+        ):  # old type1
+            self._data.phonon_dataset = self._get_dataset(self._data.phonon_supercell)
+        else:
+            self._data.phonon_dataset = self._get_dataset(
+                self._data.phonon_supercell, key_prefix="phonon_"
+            )
 
     def _parse_fc3_dataset(self):
         """Parse force dataset for fc3.
diff --git a/phono3py/version.py b/phono3py/version.py
index 00dbceca..d279b99a 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.3.3"
+__version__ = "3.3.4"