From 646a78b450eb30820c7c755816799e3adef6410b Mon Sep 17 00:00:00 2001 From: GiulioB Date: Wed, 18 Mar 2026 15:56:11 +0100 Subject: [PATCH 1/3] added broadening options gaussian_area lorentzian_height and clear names. The old gaussian and lorentzian still work but are not type hinted to discourage their usage. --- src/scm/plams/interfaces/adfsuite/ams.py | 8 +++-- src/scm/plams/tools/postprocess_results.py | 28 ++++++++++++---- unit_tests/test_tools_postprocess_results.py | 35 ++++++++++++++++++++ 3 files changed, 62 insertions(+), 9 deletions(-) create mode 100644 unit_tests/test_tools_postprocess_results.py diff --git a/src/scm/plams/interfaces/adfsuite/ams.py b/src/scm/plams/interfaces/adfsuite/ams.py index 405209d54..fc2e7aa87 100644 --- a/src/scm/plams/interfaces/adfsuite/ams.py +++ b/src/scm/plams/interfaces/adfsuite/ams.py @@ -1072,7 +1072,9 @@ 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 +1262,9 @@ 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, diff --git a/src/scm/plams/tools/postprocess_results.py b/src/scm/plams/tools/postprocess_results.py index 76c6503ad..4da3e3f04 100644 --- a/src/scm/plams/tools/postprocess_results.py +++ b/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) diff --git a/unit_tests/test_tools_postprocess_results.py b/unit_tests/test_tools_postprocess_results.py new file mode 100644 index 000000000..becfc7aef --- /dev/null +++ b/unit_tests/test_tools_postprocess_results.py @@ -0,0 +1,35 @@ +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, + ) + + integrated_area = np.trapezoid(y_result, x_result) + assert integrated_area <= np.sum(areas) From 8ca1daf6d67d86f01899b511496099fccd099742 Mon Sep 17 00:00:00 2001 From: GiulioB Date: Wed, 18 Mar 2026 16:12:14 +0100 Subject: [PATCH 2/3] mypy ams.py fixes --- src/scm/plams/interfaces/adfsuite/ams.py | 31 ++++++++++++++++++++---- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/src/scm/plams/interfaces/adfsuite/ams.py b/src/scm/plams/interfaces/adfsuite/ams.py index fc2e7aa87..d7af8284d 100644 --- a/src/scm/plams/interfaces/adfsuite/ams.py +++ b/src/scm/plams/interfaces/adfsuite/ams.py @@ -1073,7 +1073,10 @@ def get_frequency_spectrum( self, engine: Optional[str] = None, broadening_type: Literal[ - "gaussian_height", "gaussian_area", "lorentzian_height", "lorentzian_area" + "gaussian_height", + "gaussian_area", + "lorentzian_height", + "lorentzian_area", ] = "gaussian_height", broadening_width: int = 40, min_x: int = 0, @@ -1263,7 +1266,10 @@ def _get_ir_vcd_raman_spectrum( engine: Optional[str] = None, spectrum_type: Literal["ir", "raman", "vcd"] = "ir", broadening_type: Literal[ - "gaussian_height", "gaussian_area", "lorentzian_height", "lorentzian_area" + "gaussian_height", + "gaussian_area", + "lorentzian_height", + "lorentzian_area", ] = "gaussian_height", broadening_width: int = 40, min_x: int = 0, @@ -1305,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, @@ -1332,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, @@ -1359,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, From ae94ada08444234e89a31d680dff8ab76ab4a885 Mon Sep 17 00:00:00 2001 From: GiulioB Date: Wed, 18 Mar 2026 16:19:29 +0100 Subject: [PATCH 3/3] np.trapezoid is not present for python 3.8 it was called trapz --- unit_tests/test_tools_postprocess_results.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/unit_tests/test_tools_postprocess_results.py b/unit_tests/test_tools_postprocess_results.py index becfc7aef..3d377c4ba 100644 --- a/unit_tests/test_tools_postprocess_results.py +++ b/unit_tests/test_tools_postprocess_results.py @@ -30,6 +30,10 @@ def test_broaden_results_area_methods_preserve_total_area(broadening_type): broadening_width=1.0, broadening_type=broadening_type, ) - - integrated_area = np.trapezoid(y_result, x_result) + 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)