From b29a0a60198e268dac38489e4b97a40c8ba25c66 Mon Sep 17 00:00:00 2001
From: Janosh Riebesell <janosh.riebesell@gmail.com>
Date: Mon, 18 Dec 2023 16:58:32 -0800
Subject: [PATCH] Add `PhononDos.get_last_peak()` (#3525)

* add PhononDos.get_last_peak()

* add TestPhononDos.test_get_last_peak
---
 pymatgen/phonon/dos.py   | 41 ++++++++++++++++++++++++++++++++++++++++
 tests/phonon/test_dos.py |  9 +++++++++
 2 files changed, 50 insertions(+)

diff --git a/pymatgen/phonon/dos.py b/pymatgen/phonon/dos.py
index dbe5d579eea..0c06a35393e 100644
--- a/pymatgen/phonon/dos.py
+++ b/pymatgen/phonon/dos.py
@@ -354,6 +354,47 @@ def mae(self, other: PhononDos, two_sided: bool = True) -> float:
 
         return self_mae
 
+    def get_last_peak(self, threshold: float = 0.1) -> float:
+        """Find the last peak in the phonon DOS defined as the highest frequency with a DOS
+        value at least threshold * height of the overall highest DOS peak.
+        A peak is any local maximum of the DOS as a function of frequency.
+        Use dos.get_interpolated_value(peak_freq) to get density at peak_freq.
+
+        TODO method added by @janosh on 2023-12-18. seems to work well in most cases but
+        was not extensively tested. PRs with improvements welcome!
+
+        Args:
+            threshold (float, optional): Minimum ratio of the height of the last peak
+                to the height of the highest peak. Defaults to 0.1 = 10%. In case no peaks
+                are high enough to match, the threshold is reset to half the height of the
+                second-highest peak.
+
+        Returns:
+            float: last DOS peak frequency (in THz)
+        """
+        first_deriv = np.gradient(self.densities, self.frequencies)
+        second_deriv = np.gradient(first_deriv, self.frequencies)
+
+        maxima = (  # maxima indices of the first DOS derivative w.r.t. frequency
+            (first_deriv[:-1] > 0) & (first_deriv[1:] < 0) & (second_deriv[:-1] < 0)
+        )
+        # get mean of the two nearest frequencies around the maximum as better estimate
+        maxima_freqs = (self.frequencies[:-1][maxima] + self.frequencies[1:][maxima]) / 2
+
+        # filter maxima based on the threshold
+        max_dos = max(self.densities)
+        threshold = threshold * max_dos
+        filtered_maxima_freqs = maxima_freqs[self.densities[:-1][maxima] >= threshold]
+
+        if len(filtered_maxima_freqs) == 0:
+            # if no maxima reach the threshold (i.e. 1 super high peak and all other peaks
+            # tiny), use half the height of second highest peak as threshold
+            second_highest_peak = sorted(self.densities)[-2]
+            threshold = second_highest_peak / 2
+            filtered_maxima_freqs = maxima_freqs[self.densities[:-1][maxima] >= threshold]
+
+        return max(filtered_maxima_freqs)
+
 
 class CompletePhononDos(PhononDos):
     """This wrapper class defines a total dos, and also provides a list of PDos.
diff --git a/tests/phonon/test_dos.py b/tests/phonon/test_dos.py
index 757d8ed78e0..1d57946d541 100644
--- a/tests/phonon/test_dos.py
+++ b/tests/phonon/test_dos.py
@@ -102,6 +102,15 @@ def test_mae(self):
         assert self.dos.mae(dos2 - 1, two_sided=False) == pytest.approx(1.00000000000031)
         assert self.dos.mae(2 * dos2, two_sided=False) == pytest.approx(0.786546967)
 
+    def test_get_last_peak(self):
+        peak_freq = self.dos.get_last_peak()
+        assert peak_freq == approx(5.9909763191)
+        assert self.dos.get_interpolated_value(peak_freq) == approx(1.1700016497)
+
+        # try different threshold
+        peak_freq = self.dos.get_last_peak(threshold=0.5)
+        assert peak_freq == approx(4.9662820761)
+
 
 class TestCompletePhononDos(PymatgenTest):
     def setUp(self):