Degenerate heat capacity matrix

phono3py/phonon/func.py

@@ -65,8 +65,8 @@ def bose_einstein(x: NDArray[np.double], T: float) -> NDArray[np.double]:
         Temperature in K
 
     """
-    return 1.0 / (
-        np.exp(get_physical_units().THzToEv * x / (get_physical_units().KB * T)) - 1
+    return 1.0 / np.expm1(
+        get_physical_units().THzToEv * x / (get_physical_units().KB * T)
     )
 
 

phono3py/phonon/heat_capacity_matrix.py

@@ -41,17 +41,40 @@
 from phonopy.phonon.thermal_properties import mode_cv
 from phonopy.physical_units import get_physical_units
 
+# Threshold on x = E / (kB T) below which two modes are treated as
+# degenerate and (n_j - n_j') / (x_j - x_j') is replaced by its analytic
+# limit instead of the quotient.
+_DEGENERACY_TOLERANCE = 1e-6
 
-def _bose_einstein(
-    freqs: NDArray[np.double], temps: NDArray[np.double]
+
+def _bose_einstein_difference_ratio(
+    x: NDArray[np.double], n: NDArray[np.double]
 ) -> NDArray[np.double]:
-    """Bose-Einstein distribution.
+    r"""Return (n_j - n_j') / (x_j - x_j') with x = E / (k_B T).
+
+    Evaluated as a quotient, this is 0/0 for degenerate modes and loses all
+    significant digits when two frequencies differ only by round-off, which
+    happens at every q-point where bands are degenerate. The exact identity
+
+    n_{\mathbf{q}j} - n_{\mathbf{q}j'} =
+    -n_{\mathbf{q}j} (n_{\mathbf{q}j'} + 1)
+    \left( e^{x_{\mathbf{q}j} - x_{\mathbf{q}j'}} - 1 \right)
 
-    Re-implemented here because of different physical units.
+    removes the cancellation and gives the finite limit
+    -n_{\mathbf{q}j} (n_{\mathbf{q}j} + 1) at x_{\mathbf{q}j} = x_{\mathbf{q}j'}.
+    It is used for the nearly degenerate pairs, for which the exponential is
+    safely expanded to first order.
 
     """
-    x = np.divide.outer(freqs, get_physical_units().KB * temps).T
-    return 1.0 / (np.exp(x) - 1)
+    d = x[:, :, None] - x[:, None, :]
+    n_j = n[:, :, None]
+    n_jp = n[:, None, :]
+    is_deg = np.abs(d) < _DEGENERACY_TOLERANCE
+    return np.where(
+        is_deg,
+        -n_j * (n_jp + 1) * (1 + d / 2),
+        (n_j - n_jp) / np.where(is_deg, 1, d),
+    )
 
 
 def _mode_cv_matrix(
@@ -85,10 +108,11 @@ def _mode_cv_matrix(
         shape=(num_temps, num_band, num_band), dtype='double', order='C'.
 
     """
-    n = _bose_einstein(freqs, temps)
-    f_sub = np.subtract.outer(freqs, freqs)
-    n_sub = n[:, :, None] - n[:, None, :]
-    cvm = -prefactor[None, :, :] / temps[:, None, None] * n_sub / f_sub[None, :, :]
+    kb_temps = get_physical_units().KB * temps
+    x = np.divide.outer(freqs, kb_temps).T
+    n = 1.0 / np.expm1(x)
+    ratio = _bose_einstein_difference_ratio(x, n)
+    cvm = -prefactor[None, :, :] * ratio / (temps * kb_temps)[:, None, None]
     return np.ascontiguousarray(cvm)
 
 

phono3py/phonon4/api_phono4py.py

@@ -32,6 +32,7 @@
 import numpy as np
 from numpy.typing import NDArray
 from phonopy.api_phonopy import set_data_to_phonopy_yaml
+from phonopy.harmonic.dynamical_matrix import NacParams
 from phonopy.structure.atoms import PhonopyAtoms
 from phonopy.structure.cells import Primitive
 from phonopy.structure.symmetry import Symmetry
@@ -307,6 +308,7 @@ def run_frequency_shift(
         grid_points: NDArray[np.int64] | list[int] | None = None,
         temperatures: NDArray[np.double] | list[float] | None = None,
         band_indices: NDArray[np.int64] | list[int] | None = None,
+        nac_params: NacParams | None = None,
         cutoff_frequency: float = 1e-4,
         frequency_factor_to_THz: float | None = None,
     ) -> NDArray[np.double]:
@@ -325,12 +327,18 @@ def run_frequency_shift(
             Harmonic force constants defined on :attr:`supercell` (the same
             supercell fc4 is solved in), for the dynamical matrix.
         grid_points : array_like, optional
-            Grid-point indices into the no-symmetry ``np.ndindex(*mesh)``
-            addresses. Default is ``None`` (all grid points).
+            BZ-grid point indices, i.e. indices into
+            :attr:`FrequencyShift.grid_address`. Default is ``None``, which
+            uses :attr:`FrequencyShift.mesh_grid_points`, one point per mesh
+            point.
         temperatures : array_like, optional
             Temperatures in K. Default is ``[0.0]``.
         band_indices : array_like, optional
             Band indices to compute. Default is all bands.
+        nac_params : NacParams, optional
+            Non-analytical term correction parameters. Default is None, i.e.
+            no correction. The correction is dropped at Gamma, where it is
+            direction dependent and the grid gives no direction to take.
         cutoff_frequency : float, optional
             Frequencies at or below this (THz) are skipped. Default ``1e-4``.
         frequency_factor_to_THz : float, optional
@@ -358,13 +366,14 @@ def run_frequency_shift(
             mesh,
             temperatures=temperatures,
             band_indices=band_indices,
+            nac_params=nac_params,
             is_compact_fc4=is_compact_fc4,
             frequency_factor_to_THz=frequency_factor_to_THz,
             cutoff_frequency=cutoff_frequency,
             lang=self._lang,
         )
         if grid_points is None:
-            grid_points = range(int(np.prod(mesh)))
+            grid_points = self._frequency_shift.mesh_grid_points
         return np.array([self._frequency_shift.run(int(gp)) for gp in grid_points])
 
     def save(

phono3py/phonon4/frequency_shift.py

@@ -22,30 +22,18 @@
 
 import numpy as np
 from numpy.typing import NDArray
-from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
+from phonopy.harmonic.dynamical_matrix import NacParams, get_dynamical_matrix
+from phonopy.phonon.grid import BZGrid, get_qpoints_from_bz_grid_points
 from phonopy.physical_units import get_physical_units
 from phonopy.structure.atoms import PhonopyAtoms
 from phonopy.structure.cells import Primitive
 
 from phono3py._lang import resolve_lang
+from phono3py.phonon.func import bose_einstein
+from phono3py.phonon.solver import run_phonon_solver_c, run_phonon_solver_rust
 from phono3py.phonon4.real_to_reciprocal import RealToReciprocalFc4
 
 
-def _bose_einstein(frequencies: NDArray[np.double], t: float) -> NDArray[np.double]:
-    """Return Bose-Einstein occupations for frequencies (THz) at temperature t (K).
-
-    Frequencies at or below zero get zero occupation.
-    """
-    units = get_physical_units()
-    occ = np.zeros_like(frequencies)
-    if t <= 0:
-        return occ
-    mask = frequencies > 0
-    x = units.THzToEv * frequencies[mask] / (units.KB * t)
-    occ[mask] = 1.0 / (np.exp(x) - 1.0)
-    return occ
-
-
 class ReciprocalToNormalFc4:
     """Contract a reciprocal-space fc4 with eigenvectors to normal coordinates."""
 
@@ -141,6 +129,7 @@ def __init__(
         mesh: NDArray[np.int64],
         temperatures: NDArray[np.double] | None = None,
         band_indices: NDArray[np.int64] | None = None,
+        nac_params: NacParams | None = None,
         is_compact_fc4: bool = False,
         frequency_factor_to_THz: float | None = None,
         cutoff_frequency: float = 1e-4,
@@ -162,6 +151,10 @@ def __init__(
             Temperatures in K. Default is ``[0.0]``.
         band_indices : array_like, optional
             Band indices to compute. Default is all bands.
+        nac_params : NacParams, optional
+            Non-analytical term correction parameters. Default is None, i.e.
+            no correction. The correction is dropped at Gamma, where it is
+            direction dependent and the grid gives no direction to take.
         is_compact_fc4 : bool, optional
             Whether ``fc4`` is compact. Default is False.
         frequency_factor_to_THz : float, optional
@@ -188,13 +181,13 @@ def __init__(
         self._cutoff_frequency = cutoff_frequency
         self._lang = resolve_lang(lang)
 
-        self._dm = get_dynamical_matrix(fc2, supercell, primitive)
-        self._r2r = RealToReciprocalFc4(
-            fc4, primitive, self._mesh, is_compact_fc4, lang=self._lang
+        self._dm = get_dynamical_matrix(
+            fc2, supercell, primitive, nac_params=nac_params
         )
+        self._r2r = RealToReciprocalFc4(fc4, primitive, is_compact_fc4, lang=self._lang)
 
-        self._grid_address = np.array(
-            list(np.ndindex(*self._mesh.tolist())), dtype="int64"
+        self._bz_grid = BZGrid(
+            self._mesh, lattice=primitive.cell, store_dense_gp_map=True
         )
         self._frequencies, self._eigenvectors = self._solve_phonons()
         self._r2n = ReciprocalToNormalFc4(
@@ -215,57 +208,85 @@ def __init__(
             * units.EV
             / (2 * np.pi * units.THz)
             / 8
-            / np.prod(self._mesh)
+            / np.prod(self._bz_grid.D_diag)
         )
 
     def _solve_phonons(self) -> tuple[NDArray[np.double], NDArray[np.complex128]]:
-        num_grid = len(self._grid_address)
+        num_bzgp = len(self._bz_grid.addresses)
         num_band = len(self._primitive) * 3
-        frequencies = np.zeros((num_grid, num_band), dtype="double")
-        eigenvectors = np.zeros((num_grid, num_band, num_band), dtype="complex128")
-        for gi, address in enumerate(self._grid_address):
-            q = address / self._mesh
-            self._dm.run(q)
-            eigvals, eigenvectors[gi] = np.linalg.eigh(self._dm.dynamical_matrix)
-            frequencies[gi] = np.sqrt(np.abs(eigvals)) * np.sign(eigvals) * self._factor
+        frequencies = np.zeros((num_bzgp, num_band), dtype="double")
+        eigenvectors = np.zeros((num_bzgp, num_band, num_band), dtype="complex128")
+        solver = run_phonon_solver_rust if self._lang == "Rust" else run_phonon_solver_c
+        solver(
+            self._dm,
+            frequencies,
+            eigenvectors,
+            np.zeros(num_bzgp, dtype="byte"),  # phonon_done
+            np.arange(num_bzgp, dtype="int64"),
+            self._bz_grid.addresses,
+            self._bz_grid.QDinv,
+            self._factor,
+            None,  # No q-direction, so the NAC is dropped at Gamma.
+            "L",
+        )
         return frequencies, eigenvectors
 
+    @property
+    def bz_grid(self) -> BZGrid:
+        """Return the BZ-grid the phonons and the mesh sum are defined on."""
+        return self._bz_grid
+
     @property
     def grid_address(self) -> NDArray[np.int64]:
-        """Return the (no-symmetry) grid addresses."""
-        return self._grid_address
+        """Return the BZ-grid addresses."""
+        return self._bz_grid.addresses
+
+    @property
+    def mesh_grid_points(self) -> NDArray[np.int64]:
+        """Return one BZ-grid point per mesh point.
+
+        The BZ-grid holds every translationally equivalent address on the BZ
+        surface, so it has more points than the mesh. These are the unique
+        representatives the weight-one mesh sum runs over.
+
+        """
+        return self._bz_grid.grg2bzg
 
     @property
     def frequencies(self) -> NDArray[np.double]:
-        """Return phonon frequencies on the grid (THz)."""
+        """Return phonon frequencies on the BZ-grid (THz)."""
         return self._frequencies
 
     def run(self, grid_point: int) -> NDArray[np.double]:
-        """Return frequency shifts at a grid point, shape (n_temperatures, n_bands)."""
-        address0 = self._grid_address[grid_point]
-        num_grid = len(self._grid_address)
+        """Return frequency shifts at a BZ-grid point.
+
+        Returns shape (n_temperatures, n_bands).
+        """
+        q0 = get_qpoints_from_bz_grid_points(grid_point, self._bz_grid)
+        gps1 = self.mesh_grid_points
+        qpoints1 = get_qpoints_from_bz_grid_points(gps1, self._bz_grid)
         num_band = len(self._primitive) * 3
 
         # fc4_normal[gp1, band_j, band'] for the requested bands.
         fc4_normal = np.zeros(
-            (num_grid, len(self._band_indices), num_band), dtype="complex128"
+            (len(gps1), len(self._band_indices), num_band), dtype="complex128"
         )
-        for gi, address1 in enumerate(self._grid_address):
-            quartet = np.array([-address0, address0, address1, -address1])
-            fc4_reciprocal = self._r2r.run(quartet)
+        for i, (gp1, q1) in enumerate(zip(gps1, qpoints1, strict=True)):
+            fc4_reciprocal = self._r2r.run(np.array([-q0, q0, q1, -q1]))
             for j, band_index in enumerate(self._band_indices):
-                fc4_normal[gi, j] = self._r2n.run(
-                    fc4_reciprocal, grid_point, int(band_index), gi
+                fc4_normal[i, j] = self._r2n.run(
+                    fc4_reciprocal, grid_point, int(band_index), int(gp1)
                 )
 
+        freqs1 = self._frequencies[gps1]
         shifts = np.zeros((len(self._temperatures), len(self._band_indices)))
         for i_t, temperature in enumerate(self._temperatures):
-            for j in range(len(self._band_indices)):
-                total = 0.0 + 0.0j
-                for gi in range(num_grid):
-                    occ = _bose_einstein(self._frequencies[gi], temperature)
-                    total += (
-                        fc4_normal[gi, j] * self._unit_conversion * (2 * occ + 1)
-                    ).sum()
-                shifts[i_t, j] = total.real
+            occs = np.zeros_like(freqs1)
+            if temperature > 0:
+                mask = freqs1 > 0
+                occs[mask] = bose_einstein(freqs1[mask], temperature)
+            shifts[i_t] = (
+                self._unit_conversion
+                * np.einsum("gjb,gb->j", fc4_normal, 2 * occs + 1).real
+            )
         return shifts

phono3py/phonon4/real_to_reciprocal.py

@@ -29,7 +29,6 @@ def __init__(
         self,
         fc4: NDArray[np.double],
         primitive: Primitive,
-        mesh: NDArray[np.int64],
         is_compact_fc: bool = False,
         lang: Literal["C", "Rust"] = "Rust",
     ) -> None:
@@ -42,8 +41,6 @@ def __init__(
             ``(n_patom, N, N, N, 3, 3, 3, 3)`` layout.
         primitive : Primitive
             Primitive cell (provides p2s/s2p maps and smallest vectors).
-        mesh : array_like
-            Reciprocal sampling mesh, shape ``(3,)``.
         is_compact_fc : bool, optional
             Whether ``fc4`` is in the compact layout. Default is False.
         lang : {"C", "Rust"}, optional
@@ -53,7 +50,6 @@ def __init__(
         """
         self._fc4 = np.array(fc4, dtype="double", order="C")
         self._primitive = primitive
-        self._mesh = np.array(mesh, dtype="int64")
         self._is_compact_fc = is_compact_fc
         self._lang = resolve_lang(lang)
         self._p2s_map = np.array(primitive.p2s_map, dtype="int64")
@@ -64,16 +60,15 @@ def __init__(
         self._num_satom = len(self._s2p_map)
         self._num_patom = len(self._p2s_map)
 
-    def run(self, quartet: NDArray[np.int64]) -> NDArray[np.complex128]:
+    def run(self, q_points: NDArray[np.double]) -> NDArray[np.complex128]:
         """Return fc4 in reciprocal space at the given q-point quartet.
 
         Parameters
         ----------
-        quartet : ndarray
-            Four grid addresses (integers), shape ``(4, 3)``. The fractional
-            q-points are ``quartet / mesh``; the first q-point is the reference
-            (only the last three enter the phase factors), mirroring the 2015
-            convention.
+        q_points : ndarray
+            Four q-points in fractional coordinates without 2pi, shape
+            ``(4, 3)``. The first q-point is the reference (only the last three
+            enter the phase factors), mirroring the 2015 convention.
 
         Returns
         -------
@@ -83,7 +78,7 @@ def run(self, quartet: NDArray[np.int64]) -> NDArray[np.complex128]:
             dtype=complex128.
 
         """
-        q = np.array(quartet, dtype="double", order="C") / self._mesh
+        q = np.array(q_points, dtype="double", order="C")
         n_patom = self._num_patom
         fc4_reciprocal = np.zeros((n_patom,) * 4 + (3,) * 4, dtype="complex128")
         if self._lang == "Rust":