Add qpoints phonons mode as available in phonopy

Author: dmt4

docs/source/user_guide/command_line.rst

@@ -623,6 +623,13 @@ with 101 sampling points for each path segment.
    :height: 700px
    :align: center
 
+List of q-points
+++++++++++++++++
+
+Band mode outputs results in q lines or segments. Phonons data on a list of arbitrary points can be obtained with the ``--qpoints`` option. This corresponds to the ``QPOINTS`` tag in phonopy.
+
+Input q-points have to be supplied in file ``QPOINTS``, formatted as prescribed by `phonopy <https://phonopy.github.io/phonopy/input-files.html#qpoints-file>`_.
+
 
 Nudged Elastic Band
 -------------------

janus_core/calculations/phonons.py

@@ -8,7 +8,7 @@
 from ase import Atoms
 from numpy import ndarray
 import phonopy
-from phonopy.file_IO import write_force_constants_to_hdf5
+from phonopy.file_IO import parse_QPOINTS, write_force_constants_to_hdf5
 from phonopy.phonon.band_structure import (
     get_band_qpoints_and_path_connections,
     get_band_qpoints_by_seekpath,
@@ -342,11 +342,13 @@ def __init__(
         self.phonopy_file = self._build_filename("phonopy.yml")
         self.force_consts_file = self._build_filename("force_constants.hdf5")
 
-        filename = "bands" + (".hdf5" if hdf5 else ".yml")
+        suffix = ".hdf5" if hdf5 else ".yml"
+        filename = "bands" + suffix
         if not self.qpoint_file:
             filename = f"auto_{filename}"
         self.bands_file = self._build_filename(filename)
         self.bands_plot_file = self._build_filename("bands.svg")
+        self.qpoints_file = self._build_filename("qpoints" + suffix)
         self.dos_file = self._build_filename("dos.dat")
         self.dos_plot_file = self._build_filename("dos.svg")
         self.bands_dos_plot_file = self._build_filename("bs-dos.svg")
@@ -438,6 +440,11 @@ def output_files(self) -> None:
                 if self.write_results and "bands" in self.calcs
                 else None
             ),
+            "qpoints": (
+                self.qpoints_file
+                if self.write_results and "qpoints" in self.calcs
+                else None
+            ),
             "bands_plot": self.bands_plot_file if self.plot_to_file else None,
             "dos": (
                 self.dos_file if self.write_results and "dos" in self.calcs else None
@@ -712,6 +719,86 @@ def write_bands(
             build_file_dir(self.bands_plot_file)
             bplt.savefig(self.bands_plot_file)
 
+    def calc_qpoints(self, write_qpoints: bool | None = None, **kwargs) -> None:
+        """
+        Calculate phonons at qpoints supplied by file QPOINTS, analoguous to phonopy.
+
+        Parameters
+        ----------
+        write_qpoints
+            Whether to write out results to file. Default is self.write_results.
+        **kwargs
+            Additional keyword arguments to pass to `write_bands`.
+        """
+        if write_qpoints is None:
+            write_qpoints = self.write_results
+
+        # Calculate phonons if not already in results
+        if "phonon" not in self.results:
+            # Use general (self.write_results) setting for writing force constants
+            self.calc_force_constants(write_force_consts=self.write_results)
+
+        if write_qpoints:
+            self.write_qpoints(**kwargs)
+
+    def write_qpoints(
+        self,
+        *,
+        hdf5: bool | None = None,
+        qpoints_file: PathLike | None = None,
+    ) -> None:
+        """
+        Write results of qpoints mode calculations.
+
+        Parameters
+        ----------
+        hdf5
+            Whether to save the bands in an hdf5 file. Default is
+            self.hdf5.
+        qpoints_file
+            Name of yaml file to save band structure. Default is inferred from
+            `file_prefix`.
+        """
+        if "phonon" not in self.results:
+            raise ValueError(
+                "Force constants have not been calculated yet. "
+                "Please run `calc_force_constants` first"
+            )
+
+        if hdf5 is None:
+            hdf5 = self.hdf5
+
+        if qpoints_file:
+            self.qpoints_file = qpoints_file
+
+        # maybe use self.qpoint_file or allow custom input filename
+        # also allow passing a list of points programmatically
+        q_points = parse_QPOINTS()
+
+        fonons = self.results["phonon"]
+
+        fonons.run_qpoints(
+            q_points,
+            with_eigenvectors=self.write_full,
+            with_group_velocities=self.write_full,
+            with_dynamical_matrices=self.write_full,
+            # nac_q_direction = self.nac_q_direction,
+        )
+
+        build_file_dir(self.qpoints_file)
+        if hdf5:
+            # not in phonopy yet
+            # fonons.write_hdf5_qpoints_phonon(filename=self.qpoints_file)
+
+            # until the above is implemented in phonopy
+            fonons._qpoints.write_hdf5(filename=self.qpoints_file)
+        else:
+            # not in phonopy yet
+            # fonons.write_yaml_qpoints_phonon(filename=self.qpoints_file)
+
+            # until the above is implemented in phonopy
+            fonons._qpoints.write_yaml(filename=self.qpoints_file)
+
     def calc_thermal_props(
         self,
         mesh: tuple[int, int, int] | None = None,
@@ -1014,12 +1101,16 @@ def run(self) -> None:
         if "bands" in self.calcs:
             self.calc_bands()
 
+        if "qpoints" in self.calcs:
+            self.calc_qpoints()
+
         # Calculate thermal properties if specified
         if "thermal" in self.calcs:
             self.calc_thermal_props()
 
         # Calculate DOS and PDOS if specified
         if "dos" in self.calcs:
             self.calc_dos(plot_bands="bands" in self.calcs)
+
         if "pdos" in self.calcs:
             self.calc_pdos()

janus_core/cli/phonons.py

@@ -68,6 +68,13 @@ def phonons(
             help="Whether to compute band structure.", rich_help_panel="Calculation"
         ),
     ] = False,
+    qpoints: Annotated[
+        bool,
+        Option(
+            help="Whether to compute for qpoints supplied in file QPOINTS.",
+            rich_help_panel="Calculation",
+        ),
+    ] = False,
     n_qpoints: Annotated[
         int,
         Option(
@@ -218,6 +225,8 @@ def phonons(
         Mesh for sampling. Default is (10, 10, 10).
     bands
         Whether to calculate and save the band structure. Default is False.
+    qpoints
+        Whether to compute for qpoints supplied in file QPOINTS. Default is False.
     n_qpoints
         Number of q-points to sample along generated path, including end points.
         Unused if `qpoint_file` is specified. Default is 51.
@@ -360,6 +369,8 @@ def phonons(
     calcs = []
     if bands:
         calcs.append("bands")
+    if qpoints:
+        calcs.append("qpoints")
     if thermal:
         calcs.append("thermal")
     if dos:

janus_core/helpers/janus_types.py

@@ -138,7 +138,7 @@ class Correlation(TypedDict, total=True):
 Devices = Literal["cpu", "cuda", "mps", "xpu"]
 Ensembles = Literal["nph", "npt", "nve", "nvt", "nvt-nh", "nvt-csvr", "npt-mtk"]
 Properties = Literal["energy", "stress", "forces", "hessian"]
-PhononCalcs = Literal["bands", "dos", "pdos", "thermal"]
+PhononCalcs = Literal["bands", "dos", "pdos", "qpoints", "thermal"]
 Interpolators = Literal["ase", "pymatgen"]
 
 

tests/test_phonons_cli.py

@@ -146,6 +146,53 @@ def test_bands_simple(tmp_path):
     assert phonon_summary["config"]["bands"]
 
 
+def test_qpoints_simple(tmp_path):
+    """Test qpoints mode."""
+    with chdir(tmp_path):
+        file_prefix = "NaCl"
+        qpoints_results = "NaCl-qpoints.hdf5"
+        summary_path = "NaCl-phonons-summary.yml"
+
+        with open("QPOINTS", mode="w", encoding="utf8") as file:
+            file.write("3\n0.0 0.0 0.0\n0.1 0.1 0.0\n0.2 0.2 0.2\n")
+
+        result = runner.invoke(
+            app,
+            [
+                "phonons",
+                "--struct",
+                DATA_PATH / "NaCl.cif",
+                "--arch",
+                "mace_mp",
+                "--qpoints",
+                "--write-full",
+                "--hdf5",
+                "--file-prefix",
+                file_prefix,
+            ],
+        )
+        assert result.exit_code == 0
+
+        assert qpoints_results.exists()
+        with HDF5Open(qpoints_results, "r") as h5f:
+            assert "dynamical_matrix" in h5f
+            assert "eigenvector" in h5f
+            assert "frequency" in h5f
+            assert "masses" in h5f
+            assert "qpoint" in h5f
+            assert "reciprocal_lattice" in h5f
+
+        # Read phonons summary file
+        assert summary_path.exists()
+        with open(summary_path, encoding="utf8") as file:
+            phonon_summary = yaml.safe_load(file)
+
+        assert "command" in phonon_summary
+        assert "janus phonons" in phonon_summary["command"]
+        assert "config" in phonon_summary
+        assert phonon_summary["config"]["qpoints"]
+
+
 @pytest.mark.parametrize("compression", [None, "gzip", "lzf"])
 def test_hdf5(tmp_path, compression):
     """Test saving force constants and bands to HDF5 in new directory."""