Extra options for broaden_results

src/scm/plams/interfaces/adfsuite/ams.py

@@ -1072,7 +1072,12 @@ def get_frequencies(self, unit: str = "cm^-1", engine: Optional[str] = None) ->
     def get_frequency_spectrum(
         self,
         engine: Optional[str] = None,
-        broadening_type: Literal["gaussian", "lorentzian"] = "gaussian",
+        broadening_type: Literal[
+            "gaussian_height",
+            "gaussian_area",
+            "lorentzian_height",
+            "lorentzian_area",
+        ] = "gaussian_height",
         broadening_width: int = 40,
         min_x: int = 0,
         max_x: int = 4000,
@@ -1260,7 +1265,12 @@ def _get_ir_vcd_raman_spectrum(
         self,
         engine: Optional[str] = None,
         spectrum_type: Literal["ir", "raman", "vcd"] = "ir",
-        broadening_type: Literal["gaussian", "lorentzian"] = "gaussian",
+        broadening_type: Literal[
+            "gaussian_height",
+            "gaussian_area",
+            "lorentzian_height",
+            "lorentzian_area",
+        ] = "gaussian_height",
         broadening_width: int = 40,
         min_x: int = 0,
         max_x: int = 4000,
@@ -1301,7 +1311,12 @@ def _get_ir_vcd_raman_spectrum(
     def get_ir_spectrum(
         self,
         engine: Optional[str] = None,
-        broadening_type: Literal["gaussian", "lorentzian"] = "gaussian",
+        broadening_type: Literal[
+            "gaussian_height",
+            "gaussian_area",
+            "lorentzian_height",
+            "lorentzian_area",
+        ] = "gaussian_height",
         broadening_width: int = 40,
         min_x: int = 0,
         max_x: int = 4000,
@@ -1328,7 +1343,12 @@ def get_ir_spectrum(
     def get_raman_spectrum(
         self,
         engine: Optional[str] = None,
-        broadening_type: Literal["gaussian", "lorentzian"] = "gaussian",
+        broadening_type: Literal[
+            "gaussian_height",
+            "gaussian_area",
+            "lorentzian_height",
+            "lorentzian_area",
+        ] = "gaussian_height",
         broadening_width: int = 40,
         min_x: int = 0,
         max_x: int = 4000,
@@ -1355,7 +1375,12 @@ def get_raman_spectrum(
     def get_vcd_spectrum(
         self,
         engine: Optional[str] = None,
-        broadening_type: Literal["gaussian", "lorentzian"] = "gaussian",
+        broadening_type: Literal[
+            "gaussian_height",
+            "gaussian_area",
+            "lorentzian_height",
+            "lorentzian_area",
+        ] = "gaussian_height",
         broadening_width: int = 40,
         min_x: int = 0,
         max_x: int = 4000,

src/scm/plams/tools/postprocess_results.py

@@ -4,14 +4,22 @@
 ArrayOrFloat = Union[np.ndarray, float]
 
 
-def _gaussian(x: np.ndarray, A: ArrayOrFloat, x0: ArrayOrFloat, sigma: ArrayOrFloat) -> np.ndarray:
+def _gaussian_height(x: np.ndarray, A: ArrayOrFloat, x0: ArrayOrFloat, sigma: ArrayOrFloat) -> np.ndarray:
     return A * np.exp(-((x - x0) ** 2) / (2 * sigma**2))
 
 
-def _lorentzian(x: np.ndarray, A: ArrayOrFloat, x0: ArrayOrFloat, sigma: ArrayOrFloat) -> np.ndarray:
+def _lorentzian_area(x: np.ndarray, A: ArrayOrFloat, x0: ArrayOrFloat, sigma: ArrayOrFloat) -> np.ndarray:
     return (A / np.pi) * (0.5 * sigma) / ((x - x0) ** 2 + (0.5 * sigma) ** 2)
 
 
+def _gaussian_area(x: np.ndarray, A: ArrayOrFloat, x0: ArrayOrFloat, sigma: ArrayOrFloat) -> np.ndarray:
+    return (A / (sigma * (2 * np.pi) ** 0.5)) * np.exp(-((x - x0) ** 2) / (2 * sigma**2))
+
+
+def _lorentzian_height(x: np.ndarray, A: ArrayOrFloat, x0: ArrayOrFloat, sigma: ArrayOrFloat) -> np.ndarray:
+    return (A * (0.5 * sigma) ** 2) / ((x - x0) ** 2 + (0.5 * sigma) ** 2)
+
+
 def _generate_broadening_widths(x_data: np.ndarray, broadening_width: ArrayOrFloat, centers: np.ndarray) -> np.ndarray:
     if not isinstance(broadening_width, np.ndarray):
         sigmas = x_data * 0 + broadening_width
@@ -27,7 +35,9 @@ def broaden_results(
     centers: np.ndarray,
     areas: np.ndarray,
     broadening_width: Union[float, np.ndarray] = 40,
-    broadening_type: Literal["gaussian", "lorentzian"] = "gaussian",
+    broadening_type: Literal[
+        "gaussian", "lorentzian", "gaussian_height", "gaussian_area", "lorentzian_height", "lorentzian_area"
+    ] = "gaussian_height",
     x_data: Union[np.ndarray, Tuple[float, float, float]] = (0, 4000, 0.5),
     post_process: Optional[Literal["max_to_1"]] = None,
 ) -> Tuple[np.ndarray, np.ndarray]:
@@ -39,8 +49,8 @@ def broaden_results(
     :type areas: np.ndarray
     :param broadening_width: if is np.ndarray each peak broadening is assigned otherwise apply the same value for all the peaks, defaults to 40
     :type broadening_width: Union[float, np.ndarray], optional
-    :param broadening_type: the line shape of the peak, defaults to "gaussian"
-    :type broadening_type: Literal['gaussian', 'lorentzian'], optional
+    :param broadening_type: the line shape of the peak, defaults to "gaussian_height"
+    :type broadening_type: Literal['gaussian', 'lorentzian', 'gaussian_height', 'gaussian_area', 'lorentzian_height', 'lorentzian_area'], optional
     :param x_data: it can be a np.ndarray or you can set a tuple with (min,max,step_size), defaults to (0, 4000, 0.5)
     :type x_data: Union[np.ndarray, Tuple[float, float, float]], optional
     :param post_process: if max_to_1 the resulted spectrum has the max peak height =1, defaults to None
@@ -49,8 +59,12 @@ def broaden_results(
     :rtype: Tuple[np.ndarray, np.ndarray]
     """
     Function_Broaden = {
-        "gaussian": _gaussian,
-        "lorentzian": _lorentzian,
+        "gaussian": _gaussian_height,
+        "gaussian_height": _gaussian_height,
+        "lorentzian": _lorentzian_area,
+        "lorentzian_area": _lorentzian_area,
+        "lorentzian_height": _lorentzian_height,
+        "gaussian_area": _gaussian_area,
     }
     if isinstance(x_data, tuple) and len(x_data) == 3:
         x_array = np.arange(*x_data)

unit_tests/test_tools_postprocess_results.py

@@ -0,0 +1,39 @@
+import numpy as np
+import pytest
+
+from scm.plams.tools.postprocess_results import broaden_results
+
+
+@pytest.mark.parametrize("broadening_type", ["gaussian_height", "lorentzian_height", "gaussian"])
+def test_broaden_results_height_methods_have_peak_height_at_least_max_input_area(broadening_type):
+    centers = np.array([1.0, 3.0])
+    areas = np.array([1.5, 4.0])
+
+    _, y_result = broaden_results(
+        centers=centers,
+        areas=areas,
+        broadening_width=0.01,
+        broadening_type=broadening_type,
+    )
+
+    assert np.max(y_result) >= np.max(areas)
+
+
+@pytest.mark.parametrize("broadening_type", ["gaussian_area", "lorentzian_area", "lorentzian"])
+def test_broaden_results_area_methods_preserve_total_area(broadening_type):
+    centers = np.array([10.0, 30.0])
+    areas = np.array([1.5, 4.0])
+
+    x_result, y_result = broaden_results(
+        centers=centers,
+        areas=areas,
+        broadening_width=1.0,
+        broadening_type=broadening_type,
+    )
+    if hasattr(np, "trapezoid"):
+        # avoid warnings
+        integrated_area = np.trapezoid(y_result, x_result)
+    else:
+        # retrocompatibility
+        integrated_area = np.trapz(y_result, x_result)
+    assert integrated_area <= np.sum(areas)