Parse positions from vasprun.xml and use those for checking agreement…

… with displacements

phonopy/cui/create_force_sets.py

@@ -45,6 +45,8 @@
 from phonopy.interface.calculator import get_calc_dataset, get_calc_dataset_wien2k
 from phonopy.interface.lammps import rotate_lammps_forces
 from phonopy.interface.phonopy_yaml import PhonopyYaml
+from phonopy.structure.atoms import PhonopyAtoms
+from phonopy.structure.dataset import get_displacements_and_forces
 
 
 def create_FORCE_SETS(
@@ -66,8 +68,8 @@ def create_FORCE_SETS(
     """
     if log_level > 0:
         if interface_mode:
-            print("Calculator interface: %s" % interface_mode)
-        print('Displacements were read from "%s".' % disp_filename)
+            print(f"Calculator interface: {interface_mode}")
+        print(f'Displacements were read from "{disp_filename}".')
         if disp_filename == "disp.yaml":
             print("")
             print("NOTE:")
@@ -83,8 +85,8 @@ def create_FORCE_SETS(
             print("")
         if force_sets_zero_mode:
             print(
-                "Forces in %s are subtracted from forces in all "
-                "other files." % force_filenames[0]
+                f'Forces in "{force_filenames[0]}" are subtracted from forces in all '
+                "other files."
             )
 
     if disp_filename == "disp.yaml":
@@ -137,6 +139,13 @@ def create_FORCE_SETS(
             verbose=(log_level > 0),
         )
         force_sets = calc_dataset["forces"]
+        if "points" in calc_dataset:
+            if filename := _check_agreements_of_displacements(
+                supercell, disp_dataset, calc_dataset["points"], force_filenames
+            ):
+                raise RuntimeError(
+                    f'Displacements don\'t match with atomic positions in "{filename}".'
+                )
 
     if interface_mode == "lammps":
         rotate_lammps_forces(force_sets, supercell.cell, verbose=(log_level > 0))
@@ -198,6 +207,23 @@ def check_number_of_force_files(num_displacements, force_filenames, disp_filenam
         return True
 
 
+def _check_agreements_of_displacements(
+    supercell: PhonopyAtoms,
+    dataset: dict,
+    all_points: list[np.ndarray],
+    force_filenames: list[str],
+) -> Optional[str]:
+    """Check agreements of displacements."""
+    displacements = get_displacements_and_forces(dataset)[0] @ np.linalg.inv(
+        supercell.cell
+    )
+    for disp, points, filename in zip(displacements, all_points, force_filenames):
+        diff = supercell.scaled_positions + disp - points
+        diff -= np.rint(diff)
+        if (np.linalg.norm(diff @ supercell.cell, axis=1) > 1e-5).any():
+            return filename
+
+
 def _subtract_residual_forces(force_sets):
     for i in range(1, len(force_sets)):
         force_sets[i] -= force_sets[0]

phonopy/cui/phonopy_argparse.py

@@ -908,6 +908,13 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phonopy_yaml=False):
     if load_phonopy_yaml:
         parser.add_argument("filename", nargs="*", help="phonopy.yaml like file")
     else:
-        parser.add_argument("filename", nargs="*", help="Phonopy configure file")
+        parser.add_argument(
+            "filename",
+            nargs="*",
+            help=(
+                "Phonopy configure file. However if the file is recognized as "
+                "phonopy.yaml like file, this file is read as phonopy.yaml like file."
+            ),
+        )
 
     return parser, deprecated

phonopy/interface/vasp.py

@@ -145,14 +145,15 @@ def _argsort_stable(keys):
     return counts_list, reduced_symbols, sorted_positions, perm
 
 
-def _get_forces_and_energy(
+def _get_forces_points_and_energy(
     fp: io.IOBase,
     use_expat: bool = True,
     filename: Optional[Union[str, os.PathLike]] = None,
-):
+) -> tuple[np.ndarray, Optional[np.ndarray], Optional[float]]:
     vasprun = Vasprun(fp, use_expat=use_expat)
     try:
         forces = vasprun.read_forces()
+        points = vasprun.read_points()
         energy = vasprun.read_energy()
     except (RuntimeError, ValueError, xml.parsers.expat.ExpatError) as err:
         msg = (
@@ -164,7 +165,7 @@ def _get_forces_and_energy(
         if filename is not None:
             msg = f'Could not parse "{filename}". ' + msg
         raise RuntimeError(msg) from err
-    return forces, energy
+    return forces, points, energy
 
 
 def parse_set_of_forces(
@@ -179,21 +180,25 @@ def parse_set_of_forces(
 
     is_parsed = True
     force_sets = []
+    point_sets = []
     energy_sets = []
 
     for i, fp in enumerate(forces_filenames):
         if verbose:
             print(f"{i + 1}", end=" ")
         if isinstance(fp, io.IOBase):
-            forces, energy = _get_forces_and_energy(fp, use_expat=use_expat)
+            forces, points, energy = _get_forces_points_and_energy(
+                fp, use_expat=use_expat
+            )
         else:
             myio = get_io_module_to_decompress(fp)
             with myio.open(fp, "rb") as fp:
-                forces, energy = _get_forces_and_energy(
+                forces, points, energy = _get_forces_points_and_energy(
                     fp, use_expat=use_expat, filename=fp
                 )
         force_sets.append(forces)
         energy_sets.append(energy)
+        point_sets.append(points)
 
         if not check_forces(force_sets[-1], num_atoms, fp):
             is_parsed = False
@@ -202,7 +207,11 @@ def parse_set_of_forces(
         print("")
 
     if is_parsed:
-        return {"forces": force_sets, "supercell_energies": energy_sets}
+        return {
+            "forces": force_sets,
+            "points": point_sets,
+            "supercell_energies": energy_sets,
+        }
     else:
         return {}
 
@@ -613,18 +622,20 @@ def __init__(self, fileptr, use_expat=False):
         self._use_expat = use_expat
         self._vasprun_expat = None
 
-    def read_forces(self) -> Union[np.ndarray, float]:
+    def read_forces(self) -> np.ndarray:
         """Read forces either using expat or etree."""
         if self._use_expat:
-            return self._parse_expat_vasprun_xml()
+            return self._parse_expat_vasprun_xml(target="forces")
         else:
             vasprun_etree = self._parse_etree_vasprun_xml(tag="varray")
             return self._get_forces(vasprun_etree)
 
-    def read_force_constants(self):
-        """Read force constants using etree."""
-        vasprun = self._parse_etree_vasprun_xml()
-        return self._get_force_constants(vasprun)
+    def read_points(self) -> Optional[np.ndarray]:
+        """Read forces either using expat or etree."""
+        if self._use_expat:
+            return self._parse_expat_vasprun_xml(target="points")
+        else:
+            return None
 
     def read_energy(self) -> Optional[float]:
         """Read energy using expat and etree is not supported."""
@@ -633,6 +644,11 @@ def read_energy(self) -> Optional[float]:
         else:
             return None
 
+    def read_force_constants(self):
+        """Read force constants using etree."""
+        vasprun = self._parse_etree_vasprun_xml()
+        return self._get_force_constants(vasprun)
+
     def _get_forces(self, vasprun_etree):
         """Return forces using etree.
 
@@ -745,15 +761,17 @@ def _parse_by_etree(self, fileptr, tag=None):
                 yield event, elem
 
     def _parse_expat_vasprun_xml(
-        self, target: Literal["forces", "energy"] = "forces"
+        self, target: Literal["forces", "points", "energy"] = "forces"
     ) -> Union[np.ndarray, float]:
         if self._is_version528():
             return self._parse_by_expat(VasprunWrapper(self._fileptr), target=target)
         else:
             return self._parse_by_expat(self._fileptr, target=target)
 
     def _parse_by_expat(
-        self, fileptr: io.IOBase, target: Literal["forces", "energy"] = "forces"
+        self,
+        fileptr: io.IOBase,
+        target: Literal["forces", "points", "energy"] = "forces",
     ) -> Union[np.ndarray, float]:
         if self._vasprun_expat is None:
             self._vasprun_expat = VasprunxmlExpat(fileptr)
@@ -768,6 +786,8 @@ def _parse_by_expat(
 
         if target == "forces":
             return self._vasprun_expat.forces[-1]
+        if target == "points":
+            return self._vasprun_expat.points[-1]
         if target == "energy":
             return float(self._vasprun_expat.energies[-1][1])
 

test/cui/test_phonopy_cui.py

@@ -26,6 +26,7 @@ class MockArgs:
     log_level: Optional[int] = None
     fc_symmetry: bool = True
     cell_filename: Optional[str] = None
+    create_force_sets: Optional[list[str]] = None
     is_check_symmetry: Optional[bool] = None
     is_graph_plot: Optional[bool] = None
     is_graph_save: Optional[bool] = None
@@ -78,6 +79,35 @@ def test_phonopy_disp_Cr(is_ncl: bool):
         file_path.unlink()
 
 
+def test_create_force_sets():
+    """Test phonopy --create-force-sets command."""
+    argparse_control = _get_phonopy_args(
+        cell_filename=cwd / "NaCl" / "phonopy_disp.yaml.xz",
+        create_force_sets=[
+            cwd / "NaCl" / "vasprun.xml-001.xz",
+            cwd / "NaCl" / "vasprun.xml-002.xz",
+        ],
+    )
+    with pytest.raises(SystemExit) as excinfo:
+        main(**argparse_control)
+    assert excinfo.value.code == 0
+
+    for created_filename in ("FORCE_SETS",):
+        file_path = pathlib.Path(cwd_called / created_filename)
+        assert file_path.exists()
+        file_path.unlink()
+
+    argparse_control = _get_phonopy_args(
+        cell_filename=cwd / "NaCl" / "phonopy_disp.yaml.xz",
+        create_force_sets=[
+            cwd / "NaCl" / "vasprun.xml-002.xz",
+            cwd / "NaCl" / "vasprun.xml-001.xz",
+        ],
+    )
+    with pytest.raises(RuntimeError) as excinfo:
+        main(**argparse_control)
+
+
 def test_phonopy_load():
     """Test phonopy-load command."""
     # Check sys.exit(0)
@@ -114,15 +144,17 @@ def test_phonopy_is_check_symmetry():
 
 
 def _get_phonopy_args(
-    cell_filename: str,
-    supercell_dimension: Optional[str],
-    is_displacement: Optional[bool],
-    magmoms: Optional[str],
+    cell_filename: Optional[str] = None,
+    create_force_sets: Optional[list[str]] = None,
+    supercell_dimension: Optional[str] = None,
+    is_displacement: Optional[bool] = None,
+    magmoms: Optional[str] = None,
 ):
     mockargs = MockArgs(
         filename=[],
         log_level=1,
         cell_filename=cell_filename,
+        create_force_sets=create_force_sets,
         supercell_dimension=supercell_dimension,
         is_displacement=is_displacement,
         magmoms=magmoms,

test/interface/test_vasp.py

@@ -153,6 +153,12 @@ def test_parse_set_of_forces():
     dataset = parse_FORCE_SETS(filename=filename)
     force_sets = [dataset["first_atoms"][i]["forces"] for i in (0, 1)]
     energy_ref = [-216.82820693, -216.82817843]
+    np.testing.assert_allclose(
+        calc_dataset["points"][0][0], [0.00087869, 0.0, 0.0], atol=1e-5
+    )
+    np.testing.assert_allclose(
+        calc_dataset["points"][1][32], [0.25087869, 0.25, 0.25], atol=1e-5
+    )
     np.testing.assert_allclose(force_sets, calc_dataset["forces"], atol=1e-8)
     np.testing.assert_allclose(
         energy_ref, calc_dataset["supercell_energies"], atol=1e-8