From fa82ae8a4fac35029e0708adb0e1e7cf21a67695 Mon Sep 17 00:00:00 2001
From: Atsushi Togo <atz.togo@gmail.com>
Date: Sat, 4 Jan 2025 19:02:20 +0900
Subject: [PATCH 1/4] Add PbTe example

---
 example/PbTe/POSCAR-unitcell | 25 +++++++++++
 example/PbTe/README.md       | 30 +++++++++++++
 example/PbTe/velph-tmpl.toml | 85 ++++++++++++++++++++++++++++++++++++
 3 files changed, 140 insertions(+)
 create mode 100644 example/PbTe/POSCAR-unitcell
 create mode 100644 example/PbTe/README.md
 create mode 100644 example/PbTe/velph-tmpl.toml

diff --git a/example/PbTe/POSCAR-unitcell b/example/PbTe/POSCAR-unitcell
new file mode 100644
index 0000000..882f61e
--- /dev/null
+++ b/example/PbTe/POSCAR-unitcell
@@ -0,0 +1,25 @@
+generated by phonopy
+   1.00000000000000
+     6.4405958992591392    0.0000000000000000    0.0000000000000000
+     0.0000000000000000    6.4405958992591392    0.0000000000000000
+     0.0000000000000000    0.0000000000000000    6.4405958992591392
+  Pb              Te
+     4     4
+Direct
+     0.0000000000000000    0.0000000000000000    0.0000000000000000
+     0.0000000000000000    0.5000000000000000    0.5000000000000000
+     0.5000000000000000    0.0000000000000000    0.5000000000000000
+     0.5000000000000000    0.5000000000000000    0.0000000000000000
+     0.5000000000000000    0.5000000000000000    0.5000000000000000
+     0.5000000000000000    0.0000000000000000    0.0000000000000000
+     0.0000000000000000    0.5000000000000000    0.0000000000000000
+     0.0000000000000000    0.0000000000000000    0.5000000000000000
+
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
diff --git a/example/PbTe/README.md b/example/PbTe/README.md
new file mode 100644
index 0000000..a28bbea
--- /dev/null
+++ b/example/PbTe/README.md
@@ -0,0 +1,30 @@
+# PbTe example
+
+Conventional unit cell of rocksalt-type PbTe crystal structure is stored in
+`POSCAR-unitcell`. This structure is relaxed with 400 eV cutoff energy and
+PBEsol XC functional. More detail of the relaxation calculation is found in
+`velph-tmpl.toml`.
+
+2x2x2 supercell of the unit cell and the supercells with displacements are
+generated by
+
+```
+% phelel -d --dim 2 2 2 -c POSCAR-unitcell --pa auto --pm
+```
+
+By this command, following files are created. `phelel_disp.yaml` contains
+crystal structure information and `SPOSCAR`, `POSCAR-001`, `POSCAR-002`,
+`POSCAR-003`, `POSCAR-004` files are the perfect supercell and supercells with
+four different displacements, respectively.
+
+Let's assume that these supercell calculations are executed in `disp-000`,
+`disp-001`, `disp-002`, `disp-003`, and `disp-004`. After the VASP calculations
+(`LWAP=.TRUE.` or `ELPH_PREPARE=.TRUE.` in `INCAR` is necessary), phelel
+computes derivatives of local potential and PAW strengths with respect to
+displacements by
+
+```
+% phelel --cd disp-{000..004}
+```
+
+The result is stored in `phelel_params.hdf5`.
diff --git a/example/PbTe/velph-tmpl.toml b/example/PbTe/velph-tmpl.toml
new file mode 100644
index 0000000..e45eee5
--- /dev/null
+++ b/example/PbTe/velph-tmpl.toml
@@ -0,0 +1,85 @@
+[init.options]
+max_num_atoms = 600
+phono3py_max_num_atoms = 100
+
+[vasp.incar]
+encut = 400
+ncore = 4
+gga = "PS"
+
+[vasp.relax.kpoints]
+mesh = [8, 8, 8]
+shift = [0.5, 0.5, 0.5]
+
+[vasp.nac.kpoints]
+mesh = [16, 16, 16]
+
+[vasp.el_bands.kpoints]
+mesh = [16, 16, 16]
+[vasp.el_bands.kpoints_opt]
+line = 51
+[vasp.el_bands.kpoints_dense]
+mesh = [32, 32, 32]
+
+[vasp.phelel.incar]
+kpar = 2
+ediff = 1e-7
+ncore = 8
+[vasp.phelel.kpoints]
+mesh = [2, 2, 2]
+shift = [0.5, 0.5, 0.5]
+[vasp.phelel.scheduler]
+scheduler_name = "custom"
+scheduler_template = '''#!/bin/bash
+#QSUB2 core 192
+#QSUB2 mpi 192
+#QSUB2 smp 1
+#QSUB2 wtime 48:00:00
+#PBS -N {job_name}
+cd $PBS_O_WORKDIR
+
+{prepend_text}
+{mpirun_command} {vasp_binary}
+{append_text}
+'''
+job_name = "PbTe"
+mpirun_command = "mpijob"
+vasp_binary = "/home/togo/src/vasp-tmp1-atztogo/bin/vasp_std"
+prepend_text = '''
+. /etc/profile.d/modules.sh
+module load comp3h
+module load hdf5/1.14.3-inteloneapi22u3
+'''
+append_text = ""
+
+[vasp.phono3py.incar]
+kpar = 2
+addgrid = true
+[vasp.phono3py.kpoints]
+mesh = [4, 4, 4]
+shift = [0.5, 0.5, 0.5]
+[vasp.phono3py.scheduler]
+pe = "mpi* 48"
+
+[vasp.phono3py.phonon.incar]
+kpar = 2
+addgrid = true
+[vasp.phono3py.phonon.kpoints]
+mesh = [2, 2, 2]
+shift = [0.5, 0.5, 0.5]
+
+[vasp.transport.incar]
+kpar = 6
+[vasp.transport.scheduler]
+pe = "mpi* 144"
+
+[scheduler]
+scheduler_name = "sge"
+job_name = "PbTe"
+mpirun_command = "mpirun"
+vasp_binary = "/home/togo/phelel/vasp-tmp1-atztogo/bin/vasp_std"
+pe = "vienna 32"
+prepend_text = '''
+source /opt/intel/oneapi/setvars.sh --config="/home/togo/.oneapi-config"
+'''
+append_text = ""

From 5190f3f4c91292e4bd47572412148609f6188564 Mon Sep 17 00:00:00 2001
From: Atsushi Togo <atz.togo@gmail.com>
Date: Sat, 4 Jan 2025 19:04:35 +0900
Subject: [PATCH 2/4] Clear separation of phonopy instances for phelel and for
 extra phonon

---
 src/phelel/api_phelel.py                 | 226 ++++++++++++++---------
 src/phelel/base/Dij_qij.py               |   8 +-
 src/phelel/base/local_potential.py       |  44 ++---
 src/phelel/cui/load.py                   |   4 +-
 src/phelel/cui/phelel_script.py          |   4 +-
 src/phelel/file_IO.py                    |  38 ++--
 src/phelel/interface/phelel_yaml.py      |   3 +-
 src/phelel/interface/vasp/derivatives.py |  42 +++--
 test/cui/test_phelel_cui.py              |   6 +-
 test/interface/test_phelel_yaml.py       |   1 -
 test/interface/vasp/test_derivatives.py  |   3 +-
 test/phelel_disp_NaCl111.yaml            |  56 ------
 12 files changed, 218 insertions(+), 217 deletions(-)

diff --git a/src/phelel/api_phelel.py b/src/phelel/api_phelel.py
index 35b935d..ccd9d30 100644
--- a/src/phelel/api_phelel.py
+++ b/src/phelel/api_phelel.py
@@ -50,17 +50,17 @@ class Phelel:
         Phonopy class instance.
     supercell_matrix : ndarray
         Supercell matrix relative to unit cell.
-        dtype='int_', shape=(3,3)
+        dtype='long', shape=(3,3)
     primitive_matrix : ndarray
         Primitive matrix relative to unit cell. Default is None.
         dtype='double', shape=(3,3)
     atom_indices_in_derivatives : ndarray
         Atom indices in supercell used for computing derivatives.
-        dtype='int_', shape=(atom_indices_in_derivatives,)
+        dtype='long', shape=(atom_indices_in_derivatives,)
     fft_mesh : ndarray
         FFT mesh numbers for primitive cell used for generating local potential
         derivative interpolation grid.
-        dtype='int_', shape=(3,)
+        dtype='long', shape=(3,)
     dVdu : DLocalPotential
         A class instance to calculate and store dV/du.
     dDijdu : DDijQij
@@ -110,7 +110,7 @@ def __init__(
             phonon_supercell_matrix. Default is None.
         fft_mesh : array_like, optional
             FFT mesh numbers for primitive cell. Default is None.
-            dtype='int_', shape=(3,)
+            dtype='long', shape=(3,)
         symprec : float, optional
             Symmetry tolerance used to search crystal symmetry. Default
             is 1e-5.
@@ -135,30 +135,18 @@ def __init__(
         self._finufft_eps = finufft_eps
         self._log_level = log_level
 
-        ph = self._get_phonopy(supercell_matrix, primitive_matrix)
-        self._primitive_matrix = ph.primitive_matrix
-        self._supercell_matrix = ph.supercell_matrix
-        self._primitive = ph.primitive
-        self._supercell = ph.supercell
-        self._symmetry = ph.symmetry
-        self._dataset = ph.dataset
-        self._atom_indices_in_derivatives = self._primitive.p2s_map
-
-        if phonon_supercell_matrix is None:
-            self._phonon = ph
-        else:
+        self._phelel_phonon = self._get_phonopy(supercell_matrix, primitive_matrix)
+        self._atom_indices_in_derivatives = self._phelel_phonon.primitive.p2s_map
+
+        self._phonon: Optional[Phonopy] = None
+        if phonon_supercell_matrix is not None:
             self._phonon = self._get_phonopy(phonon_supercell_matrix, primitive_matrix)
-            assert isclose(self._primitive, self._phonon.primitive)
+            assert isclose(self._phelel_phonon.primitive, self._phonon.primitive)
 
         p2s_mat_float = np.linalg.inv(self.primitive.primitive_matrix)
-        self._p2s_matrix = np.rint(p2s_mat_float).astype("int_")
+        self._p2s_matrix = np.rint(p2s_mat_float).astype("long")
         assert (abs(self._p2s_matrix - p2s_mat_float) < 1e-5).all()
 
-        if phonon_supercell_matrix is None:
-            self._phonon_supercell_matrix = None
-        else:
-            self._phonon_supercell_matrix = self._phonon.supercell_matrix
-
         self._dVdu = None
         self._dDijdu = None
 
@@ -168,8 +156,8 @@ def __init__(
             self.fft_mesh = fft_mesh
 
         self._dDijdu = DDijQij(
-            self._supercell,
-            symmetry=self._symmetry,
+            self._phelel_phonon.supercell,
+            symmetry=self._phelel_phonon.symmetry,
             atom_indices=self._atom_indices_in_derivatives,
             verbose=self._log_level > 0,
         )
@@ -202,12 +190,12 @@ def unitcell(self) -> PhonopyAtoms:
     @property
     def supercell(self) -> Supercell:
         """Return supercell."""
-        return self._supercell
+        return self._phelel_phonon.supercell
 
     @property
     def primitive(self) -> Primitive:
         """Return primitive cell."""
-        return self._primitive
+        return self._phelel_phonon.primitive
 
     @property
     def p2s_matrix(self):
@@ -217,11 +205,11 @@ def p2s_matrix(self):
     @property
     def dataset(self) -> dict:
         """Setter and getter of potential dataset."""
-        return self._dataset
+        return self._phelel_phonon.dataset
 
     @dataset.setter
     def dataset(self, dataset: dict):
-        self._dataset = dataset
+        self._phelel_phonon.dataset = dataset
 
     @property
     def atom_indices_in_derivatives(self) -> np.ndarray:
@@ -231,17 +219,22 @@ def atom_indices_in_derivatives(self) -> np.ndarray:
     @property
     def symmetry(self) -> Symmetry:
         """Return symmetry of supercell."""
-        return self._symmetry
+        return self._phelel_phonon.symmetry
+
+    @property
+    def primitive_symmetry(self) -> Symmetry:
+        """Return symmetry of primitive cell."""
+        return self._phelel_phonon.primitive_symmetry
 
     @property
     def supercell_matrix(self) -> np.ndarray:
         """Return supercell matrix."""
-        return self._supercell_matrix
+        return self._phelel_phonon.supercell_matrix
 
     @property
     def primitive_matrix(self) -> Optional[np.ndarray]:
         """Return primitive matrix."""
-        return self._primitive_matrix
+        return self._phelel_phonon.primitive_matrix
 
     @property
     def phonon_supercell_matrix(self) -> Optional[np.ndarray]:
@@ -251,25 +244,39 @@ def phonon_supercell_matrix(self) -> Optional[np.ndarray]:
         ``Phelel`` class, this returns None.
 
         """
-        return self._phonon_supercell_matrix
+        if self._phonon is None:
+            return None
+        else:
+            return self._phonon.supercell_matrix
 
     @property
-    def phonon_supercell(self) -> Supercell:
+    def phonon_supercell(self) -> Optional[Supercell]:
         """Return phonon supercell."""
-        return self._phonon.supercell
+        if self._phonon is None:
+            return None
+        else:
+            return self._phonon.supercell
 
     @property
-    def phonon_primitive(self) -> Primitive:
+    def phonon_primitive(self) -> Optional[Primitive]:
         """Return phonon primitive cell."""
-        return self._phonon.primitive
+        if self._phonon is None:
+            return None
+        else:
+            return self._phonon.primitive
 
     @property
-    def phonon_dataset(self) -> dict:
+    def phonon_dataset(self) -> Optional[dict]:
         """Setter and getter of phonon dataset."""
-        return self._phonon.dataset
+        if self._phonon is None:
+            return None
+        else:
+            return self._phonon.dataset
 
     @phonon_dataset.setter
     def phonon_dataset(self, phonon_dataset):
+        if self._phonon is None:
+            raise RuntimeError("Phonon instance is not initialized.")
         self._phonon.dataset = phonon_dataset
 
     @property
@@ -288,20 +295,35 @@ def nac_params(self) -> Optional[dict]:
                 shape=(3, 3), dtype='double', order='C'
 
         """
-        return self._phonon._nac_params
+        return self._phelel_phonon.nac_params
 
     @nac_params.setter
     def nac_params(self, nac_params: dict):
-        self._phonon.nac_params = nac_params
+        self._phelel_phonon.nac_params = nac_params
 
     @property
     def force_constants(self) -> Optional[np.ndarray]:
         """Return force constants."""
-        return self._phonon.force_constants
+        if self._phonon is None:
+            return self._phelel_phonon.force_constants
+        else:
+            return self._phonon.force_constants
 
     @force_constants.setter
     def force_constants(self, force_constants: np.ndarray):
-        self._phonon.force_constants = force_constants
+        if self._phonon is None:
+            self._phelel_phonon.force_constants = force_constants
+        else:
+            self._phonon.force_constants = force_constants
+
+    @property
+    def forces(self) -> np.ndarray:
+        """Setter and getter of forces of supercells."""
+        return self._phelel_phonon.forces
+
+    @forces.setter
+    def forces(self, forces: Union[Sequence, np.ndarray]):
+        self._phelel_phonon.forces = np.array(forces, dtype="double", order="C")
 
     @property
     def unit_conversion_factor(self) -> float:
@@ -313,10 +335,13 @@ def unit_conversion_factor(self) -> float:
             (ordinary frequency).
 
         """
-        return self._phonon.unit_conversion_factor
+        if self._phonon is None:
+            return None
+        else:
+            return self._phonon.unit_conversion_factor
 
     @property
-    def phonon(self) -> Phonopy:
+    def phonon(self) -> Optional[Phonopy]:
         """Return Phonopy class instance."""
         return self._phonon
 
@@ -327,12 +352,12 @@ def fft_mesh(self) -> np.ndarray:
 
     @fft_mesh.setter
     def fft_mesh(self, fft_mesh: Union[Sequence, np.ndarray]):
-        self._fft_mesh = np.array(fft_mesh, dtype="int_")
+        self._fft_mesh = np.array(fft_mesh, dtype="long")
         self._dVdu = DLocalPotential(
             self._fft_mesh,
             self._p2s_matrix,
-            self._supercell,
-            symmetry=self._symmetry,
+            self._phelel_phonon.supercell,
+            symmetry=self._phelel_phonon.symmetry,
             atom_indices=self._atom_indices_in_derivatives,
             nufft=self._nufft,
             finufft_eps=self._finufft_eps,
@@ -366,8 +391,10 @@ def supercells_with_displacements(self) -> Optional[list[PhonopyAtoms]]:
             Phonopy.generate_displacements.
 
         """
-        ph = self._get_phonopy(self._supercell_matrix, self._primitive_matrix)
-        ph.dataset = self._dataset
+        ph = self._get_phonopy(
+            self._phelel_phonon.supercell_matrix, self._phelel_phonon.primitive_matrix
+        )
+        ph.dataset = self._phelel_phonon.dataset
         return ph.supercells_with_displacements
 
     @property
@@ -379,7 +406,10 @@ def phonon_supercells_with_displacements(self) -> Optional[list[PhonopyAtoms]]:
             Phonopy.generate_displacements.
 
         """
-        return self._phonon.supercells_with_displacements
+        if self._phonon is None:
+            return None
+        else:
+            return self._phonon.supercells_with_displacements
 
     def generate_displacements(
         self,
@@ -388,11 +418,13 @@ def generate_displacements(
         is_diagonal=True,
     ):
         """Generate displacement dataset."""
-        ph = self._get_phonopy(self._supercell_matrix, self._primitive_matrix)
+        ph = self._get_phonopy(
+            self._phelel_phonon.supercell_matrix, self._phelel_phonon.primitive_matrix
+        )
         ph.generate_displacements(
             distance=distance, is_plusminus=is_plusminus, is_diagonal=is_diagonal
         )
-        self._dataset = ph.dataset
+        self._phelel_phonon.dataset = ph.dataset
 
     def run_derivatives(self, phe_input: PhelelDataset):
         """Run displacement derivatives calculations from temporary raw data.
@@ -411,7 +443,7 @@ def run_derivatives(self, phe_input: PhelelDataset):
             raise RuntimeError(msg)
 
         if phe_input.dataset is not None:
-            self._dataset = phe_input.dataset
+            self._phelel_phonon.dataset = phe_input.dataset
         loc_pots = phe_input.local_potentials
         Dijs = phe_input.Dijs
         qijs = phe_input.qijs
@@ -424,17 +456,19 @@ def run_derivatives(self, phe_input: PhelelDataset):
             )
         else:
             self._prepare_phonon(
-                dataset=self._dataset,
+                dataset=self._phelel_phonon.dataset,
                 forces=phe_input.forces,
                 calculate_full_force_constants=True,
             )
-        self._dVdu.run(loc_pots[0], loc_pots[1:], self._dataset["first_atoms"])
+        self._dVdu.run(
+            loc_pots[0], loc_pots[1:], self._phelel_phonon.dataset["first_atoms"]
+        )
         self._dDijdu.run(
             Dijs[0],
             Dijs[1:],
             qijs[0],
             qijs[1:],
-            self._dataset["first_atoms"],
+            self._phelel_phonon.dataset["first_atoms"],
             phe_input.lm_channels,
         )
 
@@ -442,24 +476,30 @@ def save_hdf5(
         self, filename: Union[str, bytes, os.PathLike, io.IOBase] = "phelel_params.hdf5"
     ):
         """Write phelel_params.hdf5."""
-        write_phelel_params_hdf5(
-            dVdu=self._dVdu,
-            dDijdu=self._dDijdu,
-            supercell_matrix=self._supercell_matrix,
-            primitive_matrix=self._primitive_matrix,
-            primitive=self._primitive,
-            unitcell=self._unitcell,
-            supercell=self._supercell,
-            atom_indices_in_derivatives=self._atom_indices_in_derivatives,
-            disp_dataset=self._dataset,
-            force_constants=self._phonon.force_constants,
-            phonon_supercell_matrix=self._phonon.supercell_matrix,
-            phonon_primitive=self._phonon.primitive,
-            phonon_supercell=self._phonon.supercell,
-            nac_params=self._phonon.nac_params,
-            symmetry_dataset=self._phonon.primitive_symmetry.dataset,
-            filename=filename,
-        )
+        params = {
+            "dVdu": self._dVdu,
+            "dDijdu": self._dDijdu,
+            "supercell_matrix": self._phelel_phonon.supercell_matrix,
+            "primitive_matrix": self._phelel_phonon.primitive_matrix,
+            "primitive": self._phelel_phonon.primitive,
+            "unitcell": self._unitcell,
+            "supercell": self._phelel_phonon.supercell,
+            "atom_indices_in_derivatives": self._atom_indices_in_derivatives,
+            "disp_dataset": self._phelel_phonon.dataset,
+            "nac_params": self._phelel_phonon.nac_params,
+        }
+        if self._phonon is not None:
+            params.update(
+                {
+                    "force_constants": self._phonon.force_constants,
+                    "phonon_supercell_matrix": self._phonon.supercell_matrix,
+                    "phonon_primitive": self._phonon.primitive,
+                    "phonon_supercell": self._phonon.supercell,
+                    "symmetry_dataset": self._phonon.primitive_symmetry.dataset,
+                    "filename": filename,
+                }
+            )
+        write_phelel_params_hdf5(**params)
 
     def _prepare_phonon(
         self,
@@ -477,18 +517,30 @@ def _prepare_phonon(
         show_drift: bool = True,
     ):
         """Initialize phonon calculation."""
-        if dataset is not None:
-            self._phonon.dataset = dataset
-        if forces is not None:
-            self._phonon.forces = forces
-            self._phonon.produce_force_constants(
-                calculate_full_force_constants=calculate_full_force_constants,
-                fc_calculator=fc_calculator,
-                fc_calculator_options=fc_calculator_options,
-                show_drift=show_drift,
-            )
-        if force_constants is not None:
-            self._phonon.force_constants = force_constants
+        if self._phonon is not None:
+            if dataset is not None:
+                self._phonon.dataset = dataset
+            if forces is not None:
+                self._phonon.forces = forces
+                self._phonon.produce_force_constants(
+                    calculate_full_force_constants=calculate_full_force_constants,
+                    fc_calculator=fc_calculator,
+                    fc_calculator_options=fc_calculator_options,
+                    show_drift=show_drift,
+                )
+            elif force_constants is not None:
+                self._phonon.force_constants = force_constants
+        else:
+            if forces is not None:
+                self._phelel_phonon.forces = forces
+                self._phelel_phonon.produce_force_constants(
+                    calculate_full_force_constants=calculate_full_force_constants,
+                    fc_calculator=fc_calculator,
+                    fc_calculator_options=fc_calculator_options,
+                    show_drift=show_drift,
+                )
+            elif force_constants is not None:
+                self._phelel_phonon.force_constants = force_constants
 
     def _get_phonopy(
         self,
diff --git a/src/phelel/base/Dij_qij.py b/src/phelel/base/Dij_qij.py
index a808b8b..ae3e2a2 100644
--- a/src/phelel/base/Dij_qij.py
+++ b/src/phelel/base/Dij_qij.py
@@ -277,7 +277,7 @@ def _setup(self):
         else:
             atoms = self._atom_indices_in
         self._atom_indices = np.array(
-            [i for i in atoms if i in equiv_atoms], dtype="int_"
+            [i for i in atoms if i in equiv_atoms], dtype="long"
         )
         sitesym_selected_indices = [
             i for i, eq_atom in enumerate(equiv_atoms) if eq_atom in self._atom_indices
@@ -451,7 +451,7 @@ class DDijQij:
         np.unique(atom_indices given at __init__). If None, all atoms in
         supercell.
         shape=(len(atom_indices),)
-        dtype='int_'
+        dtype='long'
 
     """
 
@@ -482,9 +482,9 @@ def __init__(
         self._verbose = verbose
 
         if atom_indices is None:
-            self.atom_indices = np.arange(len(self._supercell), dtype="int_")
+            self.atom_indices = np.arange(len(self._supercell), dtype="long")
         else:
-            self.atom_indices = np.array(np.unique(atom_indices), dtype="int_")
+            self.atom_indices = np.array(np.unique(atom_indices), dtype="long")
         if symmetry is None:
             self.symmetry = Symmetry(supercell)
         else:
diff --git a/src/phelel/base/local_potential.py b/src/phelel/base/local_potential.py
index 3da379b..a5caf25 100644
--- a/src/phelel/base/local_potential.py
+++ b/src/phelel/base/local_potential.py
@@ -108,7 +108,7 @@ class LocalPotentialInterpolationNUFFT:
     ----------
     p2s_matrix : ndarray
         Transformation matrix to supercell from primitive cell
-        dtype='int_'
+        dtype='long'
         shape=(3,3)
     grid_points : ndarray
         Points to be used to compute <psi|dV_loc|psi'> given as
@@ -120,7 +120,7 @@ class LocalPotentialInterpolationNUFFT:
         shape=(prod(fft_mesh)*det(p2s_matrix), 3)
     lattice_points : ndarray
         Lattice points in supercell in primitive cell coordinates
-        dtype='int_'
+        dtype='long'
         shape=(det(supercell_matrxi), 3)
     dVdu : ndarray
         Displacement derivative of local potential in supercell interpolated on
@@ -130,7 +130,7 @@ class LocalPotentialInterpolationNUFFT:
     atom_indices_returned : ndarray
         Atom indices in supercell where dV will be computed. If those indices
         don't belong to symmetrically equivalent atoms to the dispalced atom.
-        dtype='int_', shape=(atom_indices,)
+        dtype='long', shape=(atom_indices,)
     dVs : list of DeltaLocalPotential
         This is given from DLocalPotential.run().
 
@@ -153,11 +153,11 @@ def __init__(
         ----------
         fft_mesh : array_like
             Mesh numbers for primitive cell
-            dtype='int_'
+            dtype='long'
             shape=(3,)
         p2s_matrix : array_like
             Supercell matrix relative to primitve cell
-            dtype='int_'
+            dtype='long'
             shape=(3,3)
         supercell : PhonopyAtoms
             Perfect supercell
@@ -179,7 +179,7 @@ def __init__(
         ##########
         # Inputs #
         ##########
-        self._fft_mesh = np.array(fft_mesh, dtype="int_")
+        self._fft_mesh = np.array(fft_mesh, dtype="long")
         self._verbose = verbose
         self._supercell = supercell
         self._symmetry = symmetry
@@ -196,7 +196,7 @@ def __init__(
         ##########
         # Public #
         ##########
-        self._p2s_matrix = np.array(p2s_matrix, dtype="int_", order="C")
+        self._p2s_matrix = np.array(p2s_matrix, dtype="long", order="C")
         (self._grid_points, self._lattice_points) = get_grid_points(
             fft_mesh, self._p2s_matrix
         )
@@ -327,7 +327,7 @@ def delta_Vs(self, delta_Vs: list[DeltaLocalPotential]):
         else:
             atoms = self._atom_indices_in
         self._atom_indices_returned = np.array(
-            [i for i in atoms if i in equiv_atoms], dtype="int_"
+            [i for i in atoms if i in equiv_atoms], dtype="long"
         )
         sitesym_selected_indices = [
             i
@@ -493,7 +493,7 @@ class DLocalPotential:
     ----------
     p2s_matrix : ndarray
         Transformation matrix to supercell from primitive cell.
-        dtype='int_', shape=(3,3)
+        dtype='long', shape=(3,3)
     grid_points : ndarray
         Points to be used to compute <psi|dV_loc|psi'> given as
         coordinates of supercell. Coordinates are given as
@@ -504,7 +504,7 @@ class DLocalPotential:
         shape=(prod(fft_mesh)*det(p2s_matrix), 3)
     lattice_points : ndarray
         Lattice points in supercell in primitive cell coordinates.
-        dtype='int_', shape=(det(supercell_matrxi), 3)
+        dtype='long', shape=(det(supercell_matrxi), 3)
     supercell : PhonopyAtoms
         Supercell.
     symmetry : Symmetry
@@ -517,10 +517,10 @@ class DLocalPotential:
         Atom indices in supercell where dV is computed. This is made as
         np.unique(atom_indices given at __init__). If None, all atoms in
         supercell.
-        shape=(len(atom_indices),), dtype='int_'
+        shape=(len(atom_indices),), dtype='long'
     fft_mesh : array_like
         Mesh numbers for primitive cell.
-        dtype='int_', shape=(3,)
+        dtype='long', shape=(3,)
 
     """
 
@@ -541,10 +541,10 @@ def __init__(
         ----------
         fft_mesh : array_like
             Mesh numbers for primitive cell.
-            dtype='int_', shape=(3,)
+            dtype='long', shape=(3,)
         p2s_matrix : ndarray or 2d list
             Supercell matrix relative to primitve cell.
-            dtype='int_', shape=(3,3)
+            dtype='long', shape=(3,3)
         supercell : PhonopyAtoms
             Supercell.
         symmetry : Symmetry, optional
@@ -569,9 +569,9 @@ def __init__(
         self._fft_mesh = fft_mesh
         self._p2s_matrix = p2s_matrix
         if atom_indices is None:
-            self._atom_indices = np.arange(len(self._supercell), dtype="int_")
+            self._atom_indices = np.arange(len(self._supercell), dtype="long")
         else:
-            self._atom_indices = np.array(np.unique(atom_indices), dtype="int_")
+            self._atom_indices = np.array(np.unique(atom_indices), dtype="long")
         if symmetry is None:
             self._symmetry = Symmetry(self._supercell)
         else:
@@ -644,7 +644,7 @@ def lattice_points(self):
     def lattice_points(self, lattice_points: np.ndarray):
         multi = determinant(self.p2s_matrix)
         if lattice_points.shape == (multi, 3):
-            self._lattice_points = np.array(lattice_points, dtype="int_", order="C")
+            self._lattice_points = np.array(lattice_points, dtype="long", order="C")
         else:
             raise RuntimeError(
                 "Array shape disagreement is found, %s!=%s."
@@ -862,11 +862,11 @@ def collect_site_symmetry_operations(
         Sets of indices of symmetry operations. Each set contains symmetry
         operations that send the displaced atom to specific one of equivalent
         atoms.
-        dtype='int_'
+        dtype='long'
         shape=(equiv_atoms, site_syms)
     equiv_atoms : ndarray
         Indices of symmetrically equivalent atoms
-        dtype='int_'
+        dtype='long'
         shape=(equiv_atoms,)
 
     """
@@ -881,7 +881,7 @@ def collect_site_symmetry_operations(
             for i, x in enumerate(map_to_indep_atoms)
             if map_to_indep_atoms[disp_atom] == x
         ],
-        dtype="int_",
+        dtype="long",
     )
     ops_sets = []
     for i_atom in equiv_atoms:
@@ -889,7 +889,7 @@ def collect_site_symmetry_operations(
         if len(isyms) > 0:
             ops_sets.append(isyms)
 
-    return np.array(ops_sets, dtype="int_", order="C"), equiv_atoms
+    return np.array(ops_sets, dtype="long", order="C"), equiv_atoms
 
 
 def get_displacements_with_rotations(
@@ -993,7 +993,7 @@ def get_grid_points(
         shape=(prod(fft_mesh)*det(p2s_matrix), 3)
     lattice_points : ndarray
         Lattice points in supercell in primitive cell coordinates
-        dtype='int_'
+        dtype='long'
         shape=(det(p2s_matrix), 3)
 
     """
diff --git a/src/phelel/cui/load.py b/src/phelel/cui/load.py
index f73c0e9..5feda0d 100644
--- a/src/phelel/cui/load.py
+++ b/src/phelel/cui/load.py
@@ -67,7 +67,7 @@ def load(
         FFT mesh numbers for primitive cell used for generating local potential
         derivative interpolation grid. This has to be set to run displacement
         derivative calculations not only dV/du, but also dDij/du. Default is None.
-        dtype='int_', shape=(3,)
+        dtype='long', shape=(3,)
     unitcell : PhonopyAtoms, optional
         Input unit cell. Default is None.
     supercell : PhonopyAtoms, optional
@@ -143,7 +143,7 @@ def load(
     if phonon_dataset:
         phelel.phonon_dataset = phonon_dataset
     elif dataset:
-        phelel.phonon_dataset = {k: v for k, v in dataset.items()}
+        phelel.dataset = {k: v for k, v in dataset.items()}
     if _nac_params:
         phelel.nac_params = _nac_params
 
diff --git a/src/phelel/cui/phelel_script.py b/src/phelel/cui/phelel_script.py
index 978203c..9958686 100644
--- a/src/phelel/cui/phelel_script.py
+++ b/src/phelel/cui/phelel_script.py
@@ -250,9 +250,7 @@ def main(**argparse_control):
                 phe_yml.read(filename)
 
             phelel.dataset = phe_yml.dataset
-            if phe_yml.phonon_dataset is None:
-                phelel.phonon_dataset = phe_yml.dataset
-            else:
+            if phe_yml.phonon_dataset is not None:
                 phelel.phonon_dataset = phe_yml.phonon_dataset
 
             if pathlib.Path("BORN").exists() or phe_yml.nac_params:
diff --git a/src/phelel/file_IO.py b/src/phelel/file_IO.py
index 3dbe457..4c28198 100644
--- a/src/phelel/file_IO.py
+++ b/src/phelel/file_IO.py
@@ -127,7 +127,7 @@ def read_phelel_params_hdf5(
                     np.linalg.inv(f["supercell_lattice"][:]) @ f["primitive_lattice"][:]
                 )
             p2s_mat_float = np.linalg.inv(pmat)
-            p2s_matrix = np.rint(p2s_mat_float).astype("int_")
+            p2s_matrix = np.rint(p2s_mat_float).astype("long")
             assert (abs(p2s_matrix - p2s_mat_float) < 1e-5).all()
 
         if log_level:
@@ -222,7 +222,7 @@ def _add_datasets(
         w.create_dataset("Rij", data=cmplx2real(Rij))
     if supercell_matrix is not None:
         w.create_dataset(
-            "supercell_matrix", data=np.array(supercell_matrix, dtype="int_", order="C")
+            "supercell_matrix", data=np.array(supercell_matrix, dtype="long", order="C")
         )
     if primitive_matrix is not None:
         w.create_dataset(
@@ -239,24 +239,24 @@ def _add_datasets(
             data=np.array(primitive.scaled_positions, dtype="double", order="C"),
         )
         w.create_dataset(
-            "primitive_numbers", data=np.array(primitive.numbers, dtype="int_")
+            "primitive_numbers", data=np.array(primitive.numbers, dtype="long")
         )
         w.create_dataset(
             "primitive_masses", data=np.array(primitive.masses, dtype="double")
         )
         p2s_vectors, p2s_multiplicities = _get_smallest_vectors(primitive)
-        w.create_dataset("p2s_map", data=np.array(primitive.p2s_map, dtype="int_"))
-        w.create_dataset("s2p_map", data=np.array(primitive.s2p_map, dtype="int_"))
+        w.create_dataset("p2s_map", data=np.array(primitive.p2s_map, dtype="long"))
+        w.create_dataset("s2p_map", data=np.array(primitive.s2p_map, dtype="long"))
         w.create_dataset("shortest_vectors", data=np.array(p2s_vectors, dtype="double"))
         w.create_dataset(
             "shortest_vector_multiplicities",
-            data=np.array(p2s_multiplicities, dtype="int_"),
+            data=np.array(p2s_multiplicities, dtype="long"),
         )
         if atom_indices_in_derivatives is not None:
             if not np.array_equal(atom_indices_in_derivatives, primitive.p2s_map):
                 w.create_dataset(
                     "atom_indices_in_derivatives",
-                    data=np.array(atom_indices_in_derivatives, dtype="int_"),
+                    data=np.array(atom_indices_in_derivatives, dtype="long"),
                 )
     if unitcell is not None:
         w.create_dataset(
@@ -268,7 +268,7 @@ def _add_datasets(
             data=np.array(unitcell.scaled_positions, dtype="double", order="C"),
         )
         w.create_dataset(
-            "unitcell_numbers", data=np.array(unitcell.numbers, dtype="int_")
+            "unitcell_numbers", data=np.array(unitcell.numbers, dtype="long")
         )
         w.create_dataset(
             "unitcell_masses", data=np.array(unitcell.masses, dtype="double")
@@ -283,7 +283,7 @@ def _add_datasets(
             data=np.array(supercell.scaled_positions, dtype="double", order="C"),
         )
         w.create_dataset(
-            "supercell_numbers", data=np.array(supercell.numbers, dtype="int_")
+            "supercell_numbers", data=np.array(supercell.numbers, dtype="long")
         )
         w.create_dataset(
             "supercell_masses", data=np.array(supercell.masses, dtype="double")
@@ -292,7 +292,7 @@ def _add_datasets(
         if "first_atoms" in disp_dataset:
             atom_indices = [d["number"] for d in disp_dataset["first_atoms"]]
             w.create_dataset(
-                "displacements_atom_indices", data=np.array(atom_indices, dtype="int_")
+                "displacements_atom_indices", data=np.array(atom_indices, dtype="long")
             )
             disps = [d["displacement"] for d in disp_dataset["first_atoms"]]
             w.create_dataset(
@@ -305,25 +305,25 @@ def _add_datasets(
     if phonon_supercell_matrix is not None:
         w.create_dataset(
             "phonon_supercell_matrix",
-            data=np.array(phonon_supercell_matrix, dtype="int_", order="C"),
+            data=np.array(phonon_supercell_matrix, dtype="long", order="C"),
         )
     if phonon_primitive is not None:
         p2s_vectors, p2s_multiplicities = _get_smallest_vectors(phonon_primitive)
         w.create_dataset(
-            "phonon_p2s_map", data=np.array(phonon_primitive.p2s_map, dtype="int_")
+            "phonon_p2s_map", data=np.array(phonon_primitive.p2s_map, dtype="long")
         )
         w.create_dataset(
-            "phonon_s2p_map", data=np.array(phonon_primitive.s2p_map, dtype="int_")
+            "phonon_s2p_map", data=np.array(phonon_primitive.s2p_map, dtype="long")
         )
         w.create_dataset(
             "phonon_shortest_vectors", data=np.array(p2s_vectors, dtype="double")
         )
         w.create_dataset(
             "phonon_shortest_vector_multiplicities",
-            data=np.array(p2s_multiplicities, dtype="int_"),
+            data=np.array(p2s_multiplicities, dtype="long"),
         )
         p2s_mat_float = np.linalg.inv(phonon_primitive.primitive_matrix)
-        p2s_matrix = np.rint(p2s_mat_float).astype("int_")
+        p2s_matrix = np.rint(p2s_mat_float).astype("long")
         assert (abs(p2s_matrix - p2s_mat_float) < 1e-5).all()
         lattice_points, _ = get_lattice_points(p2s_matrix)
         w.create_dataset("phonon_lattice_point", data=lattice_points)
@@ -338,7 +338,7 @@ def _add_datasets(
         )
         w.create_dataset(
             "phonon_supercell_numbers",
-            data=np.array(phonon_supercell.numbers, dtype="int_"),
+            data=np.array(phonon_supercell.numbers, dtype="long"),
         )
         w.create_dataset(
             "phonon_supercell_masses",
@@ -362,18 +362,18 @@ def _add_datasets(
         )
         w.create_dataset(
             "direct_rotations",
-            data=np.array(symmetry_dataset.rotations, dtype="int_", order="C"),
+            data=np.array(symmetry_dataset.rotations, dtype="long", order="C"),
         )
         sym_dataset_dict = vars(symmetry_dataset)
         if "number" in sym_dataset_dict:  # For non-magnetic case
             w.create_dataset(
-                "spacegroup_number", data=int(symmetry_dataset.number), dtype="int_"
+                "spacegroup_number", data=int(symmetry_dataset.number), dtype="long"
             )
         elif "uni_number" in sym_dataset_dict:  # For magnetic case
             w.create_dataset(
                 "magnetic_spacegroup_uni_number",
                 data=int(symmetry_dataset.uni_number),
-                dtype="int_",
+                dtype="long",
             )
 
 
diff --git a/src/phelel/interface/phelel_yaml.py b/src/phelel/interface/phelel_yaml.py
index 0409085..60633b3 100644
--- a/src/phelel/interface/phelel_yaml.py
+++ b/src/phelel/interface/phelel_yaml.py
@@ -189,7 +189,8 @@ def _nac_yaml_lines(self):
         This method override PhonopyYaml._nac_yaml_lines.
 
         """
-        assert isclose(self._data.primitive, self._data.phonon_primitive)
+        if self._data.phonon_primitive is not None:
+            assert isclose(self._data.primitive, self._data.phonon_primitive)
         return super()._nac_yaml_lines()
 
     def _displacements_yaml_lines(self, with_forces=False):
diff --git a/src/phelel/interface/vasp/derivatives.py b/src/phelel/interface/vasp/derivatives.py
index 7fe3dec..1151062 100644
--- a/src/phelel/interface/vasp/derivatives.py
+++ b/src/phelel/interface/vasp/derivatives.py
@@ -8,9 +8,11 @@
 from typing import Optional, Union
 
 import numpy as np
-from phonopy import Phonopy
 from phonopy.file_IO import get_born_parameters
 from phonopy.interface.vasp import parse_set_of_forces
+from phonopy.structure.atoms import PhonopyAtoms
+from phonopy.structure.cells import Primitive
+from phonopy.structure.symmetry import Symmetry
 
 from phelel import Phelel
 from phelel.api_phelel import PhelelDataset
@@ -58,23 +60,37 @@ def read_files(
     else:
         _dir_names = phonon_dir_names
     vasprun_filenames = _get_vasprun_filenames(_dir_names, log_level=log_level)
+
+    if phelel.phonon_supercell_matrix:
+        supercell = phelel.phonon_supercell
+    else:
+        supercell = phelel.supercell
     forces = read_forces_from_vasprunxmls(
         vasprun_filenames,
-        phelel.phonon,
+        supercell,
         subtract_rfs=subtract_rfs,
         log_level=log_level,
     )
 
-    if forces[0].shape[0] != len(phelel.phonon_supercell):
+    if forces[0].shape[0] != len(supercell):
         raise ValueError(
             "Number of ions in the phonon supercell is different from the number of "
             "atoms in the vasprun.xml file."
         )
 
-    nac_params = _read_born(phelel.phonon)
-    phelel.phonon.forces = forces
+    if phelel.phonon_supercell_matrix is None:
+        phelel.forces = forces
+    else:
+        phelel.phonon.forces = forces
+
+    nac_params = _read_born(
+        phelel.primitive, phelel.primitive_symmetry, log_level=log_level
+    )
     if nac_params:
-        phelel.phonon.nac_params = nac_params
+        phelel.nac_params = nac_params
+        if phelel.phonon_supercell_matrix is not None:
+            phelel.phonon.nac_params = nac_params
+
     return PhelelDataset(
         local_potentials=loc_pots,
         Dijs=Dijs,
@@ -138,14 +154,12 @@ def read_Rij(dir_name, inwap_per):
 
 def read_forces_from_vasprunxmls(
     vasprun_filenames: Union[list, tuple],
-    phonon: Phonopy,
+    supercell: PhonopyAtoms,
     subtract_rfs=False,
     log_level=0,
 ):
     """Read forces from vasprun.xml's and read NAC params from BORN."""
-    calc_dataset = parse_set_of_forces(
-        len(phonon.supercell), vasprun_filenames, verbose=False
-    )
+    calc_dataset = parse_set_of_forces(len(supercell), vasprun_filenames, verbose=False)
     forces = calc_dataset["forces"]
 
     if subtract_rfs:
@@ -171,12 +185,10 @@ def _get_vasprun_filenames(dir_names, log_level=0):
     return vasprun_filenames
 
 
-def _read_born(phonon: Phonopy, log_level=0):
+def _read_born(primitive: Primitive, primitive_symmetry: Symmetry, log_level: int = 0):
     if pathlib.Path("BORN").exists():
         with open("BORN", "r") as f:
-            nac_params = get_born_parameters(
-                f, phonon.primitive, phonon.primitive_symmetry
-            )
+            nac_params = get_born_parameters(f, primitive, primitive_symmetry)
             if log_level:
                 print('"BORN" was read.')
         return nac_params
@@ -188,7 +200,7 @@ def _get_datasets(phelel: Phelel) -> tuple:
     """Return inwap dataset and phonopy dataset."""
     if "first_atoms" in phelel.dataset:
         dataset = phelel.dataset
-        if "first_atoms" in phelel.phonon_dataset:
+        if phelel.phonon_supercell_matrix and "first_atoms" in phelel.phonon_dataset:
             phonon_dataset = phelel.phonon_dataset
         else:
             phonon_dataset = dataset
diff --git a/test/cui/test_phelel_cui.py b/test/cui/test_phelel_cui.py
index aafe2d8..a19b2f6 100644
--- a/test/cui/test_phelel_cui.py
+++ b/test/cui/test_phelel_cui.py
@@ -124,8 +124,6 @@ def test_phelel_script_create_displacements(is_plusminus_displacements: bool):
             "phelel_disp.yaml",
             "SPOSCAR",
             "POSCAR-002",
-            "SPOSCAR_PH",
-            "POSCAR_PH-002",
         )
         for created_filename in ("POSCAR-003", "POSCAR_PH-003"):
             file_path = pathlib.Path(cwd_called / created_filename)
@@ -135,10 +133,8 @@ def test_phelel_script_create_displacements(is_plusminus_displacements: bool):
             "phelel_disp.yaml",
             "SPOSCAR",
             "POSCAR-001",
-            "SPOSCAR_PH",
-            "POSCAR_PH-001",
         )
-        for created_filename in ("POSCAR-002", "POSCAR_PH-002"):
+        for created_filename in ("POSCAR-002",):
             file_path = pathlib.Path(cwd_called / created_filename)
             assert not file_path.exists()
 
diff --git a/test/interface/test_phelel_yaml.py b/test/interface/test_phelel_yaml.py
index cb207b9..e00f927 100644
--- a/test/interface/test_phelel_yaml.py
+++ b/test/interface/test_phelel_yaml.py
@@ -23,7 +23,6 @@ def test_PhelelYaml_get_yaml_lines():
     assert isclose(phe_yml.primitive, phe_yml_test.primitive)
     assert isclose(phe_yml.unitcell, phe_yml_test.unitcell)
     assert isclose(phe_yml.supercell, phe_yml_test.supercell)
-    assert isclose(phe_yml.phonon_supercell, phe_yml_test.phonon_supercell)
     assert phe_yml.version == phe_yml_test.version
     np.testing.assert_array_equal(
         phe_yml.supercell_matrix, phe_yml_test.supercell_matrix
diff --git a/test/interface/vasp/test_derivatives.py b/test/interface/vasp/test_derivatives.py
index a017423..7b6b3e5 100644
--- a/test/interface/vasp/test_derivatives.py
+++ b/test/interface/vasp/test_derivatives.py
@@ -23,7 +23,7 @@ def test_read_forces_from_vasprunxmls_NaCl111(phelel_empty_NaCl111: Phelel):
     ]
     forces = read_forces_from_vasprunxmls(
         filenames,
-        phe.phonon,
+        phe.supercell,
         subtract_rfs=True,
         log_level=0,
     )
@@ -103,5 +103,4 @@ def _get_phelel_C111(phelel_yaml_filename: str) -> Phelel:
         log_level=1,
     )
     phelel.dataset = phe_yml.dataset
-    phelel.phonon.dataset = phe_yml.dataset
     return phelel
diff --git a/test/phelel_disp_NaCl111.yaml b/test/phelel_disp_NaCl111.yaml
index fbcb54b..db83a2a 100644
--- a/test/phelel_disp_NaCl111.yaml
+++ b/test/phelel_disp_NaCl111.yaml
@@ -128,62 +128,6 @@ supercell:
     mass: 35.453000
     reduced_to: 5
 
-phonon_primitive_cell:
-  lattice:
-  - [     0.000000000000000,     2.845150738087836,     2.845150738087836 ] # a
-  - [     2.845150738087836,     0.000000000000000,     2.845150738087836 ] # b
-  - [     2.845150738087836,     2.845150738087836,     0.000000000000000 ] # c
-  points:
-  - symbol: Na # 1
-    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
-    mass: 22.989769
-  - symbol: Cl # 2
-    coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
-    mass: 35.453000
-  reciprocal_lattice: # without 2pi
-  - [    -0.175737613233118,     0.175737613233118,     0.175737613233118 ] # a*
-  - [     0.175737613233118,    -0.175737613233118,     0.175737613233118 ] # b*
-  - [     0.175737613233118,     0.175737613233118,    -0.175737613233118 ] # c*
-
-phonon_supercell:
-  lattice:
-  - [     5.690301476175671,     0.000000000000000,     0.000000000000000 ] # a
-  - [     0.000000000000000,     5.690301476175671,     0.000000000000000 ] # b
-  - [     0.000000000000000,     0.000000000000000,     5.690301476175671 ] # c
-  points:
-  - symbol: Na # 1
-    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
-    mass: 22.989769
-    reduced_to: 1
-  - symbol: Na # 2
-    coordinates: [  0.000000000000000,  0.500000000000000,  0.500000000000000 ]
-    mass: 22.989769
-    reduced_to: 1
-  - symbol: Na # 3
-    coordinates: [  0.500000000000000,  0.000000000000000,  0.500000000000000 ]
-    mass: 22.989769
-    reduced_to: 1
-  - symbol: Na # 4
-    coordinates: [  0.500000000000000,  0.500000000000000,  0.000000000000000 ]
-    mass: 22.989769
-    reduced_to: 1
-  - symbol: Cl # 5
-    coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
-    mass: 35.453000
-    reduced_to: 5
-  - symbol: Cl # 6
-    coordinates: [  0.500000000000000,  0.000000000000000,  0.000000000000000 ]
-    mass: 35.453000
-    reduced_to: 5
-  - symbol: Cl # 7
-    coordinates: [  0.000000000000000,  0.500000000000000,  0.000000000000000 ]
-    mass: 35.453000
-    reduced_to: 5
-  - symbol: Cl # 8
-    coordinates: [  0.000000000000000,  0.000000000000000,  0.500000000000000 ]
-    mass: 35.453000
-    reduced_to: 5
-
 displacements:
 - atom:    1
   displacement:

From 1d592ecfb862b1bf959b755a7d12da88a197213c Mon Sep 17 00:00:00 2001
From: Atsushi Togo <atz.togo@gmail.com>
Date: Sat, 4 Jan 2025 19:23:06 +0900
Subject: [PATCH 3/4] Further refactoring of Phelel class

---
 src/phelel/api_phelel.py | 64 ++++++++++++++++++----------------------
 1 file changed, 29 insertions(+), 35 deletions(-)

diff --git a/src/phelel/api_phelel.py b/src/phelel/api_phelel.py
index ccd9d30..a531146 100644
--- a/src/phelel/api_phelel.py
+++ b/src/phelel/api_phelel.py
@@ -300,6 +300,8 @@ def nac_params(self) -> Optional[dict]:
     @nac_params.setter
     def nac_params(self, nac_params: dict):
         self._phelel_phonon.nac_params = nac_params
+        if self._phonon is not None:
+            self._phonon.nac_params = nac_params
 
     @property
     def force_constants(self) -> Optional[np.ndarray]:
@@ -319,11 +321,18 @@ def force_constants(self, force_constants: np.ndarray):
     @property
     def forces(self) -> np.ndarray:
         """Setter and getter of forces of supercells."""
-        return self._phelel_phonon.forces
+        if self._phonon is None:
+            return self._phelel_phonon.forces
+        else:
+            return self._phonon.forces
 
     @forces.setter
     def forces(self, forces: Union[Sequence, np.ndarray]):
-        self._phelel_phonon.forces = np.array(forces, dtype="double", order="C")
+        _forces = np.array(forces, dtype="double", order="C")
+        if self._phonon is None:
+            self._phelel_phonon.forces = _forces
+        else:
+            self._phonon.forces = _forces
 
     @property
     def unit_conversion_factor(self) -> float:
@@ -391,11 +400,7 @@ def supercells_with_displacements(self) -> Optional[list[PhonopyAtoms]]:
             Phonopy.generate_displacements.
 
         """
-        ph = self._get_phonopy(
-            self._phelel_phonon.supercell_matrix, self._phelel_phonon.primitive_matrix
-        )
-        ph.dataset = self._phelel_phonon.dataset
-        return ph.supercells_with_displacements
+        return self._phelel_phonon.supercells_with_displacements
 
     @property
     def phonon_supercells_with_displacements(self) -> Optional[list[PhonopyAtoms]]:
@@ -418,13 +423,9 @@ def generate_displacements(
         is_diagonal=True,
     ):
         """Generate displacement dataset."""
-        ph = self._get_phonopy(
-            self._phelel_phonon.supercell_matrix, self._phelel_phonon.primitive_matrix
-        )
-        ph.generate_displacements(
+        self._phelel_phonon.generate_displacements(
             distance=distance, is_plusminus=is_plusminus, is_diagonal=is_diagonal
         )
-        self._phelel_phonon.dataset = ph.dataset
 
     def run_derivatives(self, phe_input: PhelelDataset):
         """Run displacement derivatives calculations from temporary raw data.
@@ -517,30 +518,23 @@ def _prepare_phonon(
         show_drift: bool = True,
     ):
         """Initialize phonon calculation."""
-        if self._phonon is not None:
-            if dataset is not None:
-                self._phonon.dataset = dataset
-            if forces is not None:
-                self._phonon.forces = forces
-                self._phonon.produce_force_constants(
-                    calculate_full_force_constants=calculate_full_force_constants,
-                    fc_calculator=fc_calculator,
-                    fc_calculator_options=fc_calculator_options,
-                    show_drift=show_drift,
-                )
-            elif force_constants is not None:
-                self._phonon.force_constants = force_constants
+        if self._phonon is None:
+            ph = self._phelel_phonon
         else:
-            if forces is not None:
-                self._phelel_phonon.forces = forces
-                self._phelel_phonon.produce_force_constants(
-                    calculate_full_force_constants=calculate_full_force_constants,
-                    fc_calculator=fc_calculator,
-                    fc_calculator_options=fc_calculator_options,
-                    show_drift=show_drift,
-                )
-            elif force_constants is not None:
-                self._phelel_phonon.force_constants = force_constants
+            ph = self._phonon
+        if dataset is not None:
+            ph.dataset = dataset
+        if forces is not None:
+            ph.forces = forces
+        if ph.forces is not None:
+            ph.produce_force_constants(
+                calculate_full_force_constants=calculate_full_force_constants,
+                fc_calculator=fc_calculator,
+                fc_calculator_options=fc_calculator_options,
+                show_drift=show_drift,
+            )
+        elif force_constants is not None:
+            ph.force_constants = force_constants
 
     def _get_phonopy(
         self,

From b45f30fbb58ca420b4f64d2bf439ad0838806c82 Mon Sep 17 00:00:00 2001
From: Atsushi Togo <atz.togo@gmail.com>
Date: Sat, 4 Jan 2025 19:24:33 +0900
Subject: [PATCH 4/4] Add example phelel_disp.yaml

---
 example/PbTe/phelel_disp_reference.yaml | 362 ++++++++++++++++++++++++
 1 file changed, 362 insertions(+)
 create mode 100644 example/PbTe/phelel_disp_reference.yaml

diff --git a/example/PbTe/phelel_disp_reference.yaml b/example/PbTe/phelel_disp_reference.yaml
new file mode 100644
index 0000000..b3c884d
--- /dev/null
+++ b/example/PbTe/phelel_disp_reference.yaml
@@ -0,0 +1,362 @@
+phelel:
+  version: "0.7.0"
+  calculator: vasp
+  symmetry_tolerance: 1.00000e-05
+  configuration:
+    cell_filename: "POSCAR-unitcell"
+    create_displacements: ".true."
+    pm: ".true."
+    primitive_axes: "auto"
+    dim: "2 2 2"
+
+physical_unit:
+  atomic_mass: "AMU"
+  length: "angstrom"
+
+space_group:
+  type: "Fm-3m"
+  number: 225
+  Hall_symbol: "-F 4 2 3"
+
+primitive_matrix:
+- [  0.000000000000000,  0.500000000000000,  0.500000000000000 ]
+- [  0.500000000000000,  0.000000000000000,  0.500000000000000 ]
+- [  0.500000000000000,  0.500000000000000,  0.000000000000000 ]
+
+supercell_matrix:
+- [   2,   0,   0 ]
+- [   0,   2,   0 ]
+- [   0,   0,   2 ]
+
+primitive_cell:
+  lattice:
+  - [     0.000000000000000,     3.220297949629569,     3.220297949629569 ] # a
+  - [     3.220297949629569,     0.000000000000000,     3.220297949629569 ] # b
+  - [     3.220297949629569,     3.220297949629569,     0.000000000000000 ] # c
+  points:
+  - symbol: Pb # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 207.200000
+  - symbol: Te # 2
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 127.600000
+  reciprocal_lattice: # without 2pi
+  - [    -0.155265136276448,     0.155265136276448,     0.155265136276448 ] # a*
+  - [     0.155265136276448,    -0.155265136276448,     0.155265136276448 ] # b*
+  - [     0.155265136276448,     0.155265136276448,    -0.155265136276448 ] # c*
+
+unit_cell:
+  lattice:
+  - [     6.440595899259139,     0.000000000000000,     0.000000000000000 ] # a
+  - [     0.000000000000000,     6.440595899259139,     0.000000000000000 ] # b
+  - [     0.000000000000000,     0.000000000000000,     6.440595899259139 ] # c
+  points:
+  - symbol: Pb # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 2
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 3
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 4
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Te # 5
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 5
+  - symbol: Te # 6
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 5
+  - symbol: Te # 7
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 5
+  - symbol: Te # 8
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 5
+
+supercell:
+  lattice:
+  - [    12.881191798518278,     0.000000000000000,     0.000000000000000 ] # a
+  - [     0.000000000000000,    12.881191798518278,     0.000000000000000 ] # b
+  - [     0.000000000000000,     0.000000000000000,    12.881191798518278 ] # c
+  points:
+  - symbol: Pb # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 2
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 3
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 4
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 5
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 6
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 7
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 8
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 9
+    coordinates: [  0.000000000000000,  0.250000000000000,  0.250000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 10
+    coordinates: [  0.500000000000000,  0.250000000000000,  0.250000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 11
+    coordinates: [  0.000000000000000,  0.750000000000000,  0.250000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 12
+    coordinates: [  0.500000000000000,  0.750000000000000,  0.250000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 13
+    coordinates: [  0.000000000000000,  0.250000000000000,  0.750000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 14
+    coordinates: [  0.500000000000000,  0.250000000000000,  0.750000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 15
+    coordinates: [  0.000000000000000,  0.750000000000000,  0.750000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 16
+    coordinates: [  0.500000000000000,  0.750000000000000,  0.750000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 17
+    coordinates: [  0.250000000000000,  0.000000000000000,  0.250000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 18
+    coordinates: [  0.750000000000000,  0.000000000000000,  0.250000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 19
+    coordinates: [  0.250000000000000,  0.500000000000000,  0.250000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 20
+    coordinates: [  0.750000000000000,  0.500000000000000,  0.250000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 21
+    coordinates: [  0.250000000000000,  0.000000000000000,  0.750000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 22
+    coordinates: [  0.750000000000000,  0.000000000000000,  0.750000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 23
+    coordinates: [  0.250000000000000,  0.500000000000000,  0.750000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 24
+    coordinates: [  0.750000000000000,  0.500000000000000,  0.750000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 25
+    coordinates: [  0.250000000000000,  0.250000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 26
+    coordinates: [  0.750000000000000,  0.250000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 27
+    coordinates: [  0.250000000000000,  0.750000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 28
+    coordinates: [  0.750000000000000,  0.750000000000000,  0.000000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 29
+    coordinates: [  0.250000000000000,  0.250000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 30
+    coordinates: [  0.750000000000000,  0.250000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 31
+    coordinates: [  0.250000000000000,  0.750000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Pb # 32
+    coordinates: [  0.750000000000000,  0.750000000000000,  0.500000000000000 ]
+    mass: 207.200000
+    reduced_to: 1
+  - symbol: Te # 33
+    coordinates: [  0.250000000000000,  0.250000000000000,  0.250000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 34
+    coordinates: [  0.750000000000000,  0.250000000000000,  0.250000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 35
+    coordinates: [  0.250000000000000,  0.750000000000000,  0.250000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 36
+    coordinates: [  0.750000000000000,  0.750000000000000,  0.250000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 37
+    coordinates: [  0.250000000000000,  0.250000000000000,  0.750000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 38
+    coordinates: [  0.750000000000000,  0.250000000000000,  0.750000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 39
+    coordinates: [  0.250000000000000,  0.750000000000000,  0.750000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 40
+    coordinates: [  0.750000000000000,  0.750000000000000,  0.750000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 41
+    coordinates: [  0.250000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 42
+    coordinates: [  0.750000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 43
+    coordinates: [  0.250000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 44
+    coordinates: [  0.750000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 45
+    coordinates: [  0.250000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 46
+    coordinates: [  0.750000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 47
+    coordinates: [  0.250000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 48
+    coordinates: [  0.750000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 49
+    coordinates: [  0.000000000000000,  0.250000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 50
+    coordinates: [  0.500000000000000,  0.250000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 51
+    coordinates: [  0.000000000000000,  0.750000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 52
+    coordinates: [  0.500000000000000,  0.750000000000000,  0.000000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 53
+    coordinates: [  0.000000000000000,  0.250000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 54
+    coordinates: [  0.500000000000000,  0.250000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 55
+    coordinates: [  0.000000000000000,  0.750000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 56
+    coordinates: [  0.500000000000000,  0.750000000000000,  0.500000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 57
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.250000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 58
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.250000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 59
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.250000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 60
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.250000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 61
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.750000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 62
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.750000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 63
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.750000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+  - symbol: Te # 64
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.750000000000000 ]
+    mass: 127.600000
+    reduced_to: 33
+
+displacements:
+- atom:    1
+  displacement:
+    [   0.0300000000000000,  0.0000000000000000,  0.0000000000000000 ]
+- atom:    1
+  displacement:
+    [  -0.0300000000000000,  0.0000000000000000,  0.0000000000000000 ]
+- atom:   33
+  displacement:
+    [   0.0300000000000000,  0.0000000000000000,  0.0000000000000000 ]
+- atom:   33
+  displacement:
+    [  -0.0300000000000000,  0.0000000000000000,  0.0000000000000000 ]