From 474100090e19cc55dd332c79459e5a6b476cc0f6 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Mon, 24 Aug 2026 09:20:43 +0200 Subject: [PATCH] fix: validate rdf selections, delta_r and cell volume The RDF setup accepted inputs it cannot handle: a selection matching no atoms produced an all-NaN g(r) on the general path and a bare ZeroDivisionError on the legacy path, delta_r = 0.0 aborted with an unhandled OverflowError or ZeroDivisionError, and a vacuum trajectory with an explicit radial range divided by an infinite average volume and silently wrote NaN and inf into every output column. Each case now raises a clear RDFError during setup, matching the empty selection handling of the MSD, VACF and momentum analyses. The user guide no longer claims vacuum trajectories are supported. --- PQAnalysis/analysis/rdf/rdf.py | 47 ++++++++- docs/source/userGuide/userGuide.rst | 5 +- tests/analysis/rdf/test_rdf.py | 149 +++++++++++++++++++++++++++- 3 files changed, 193 insertions(+), 8 deletions(-) diff --git a/PQAnalysis/analysis/rdf/rdf.py b/PQAnalysis/analysis/rdf/rdf.py index ff8e37d1..98f92af6 100644 --- a/PQAnalysis/analysis/rdf/rdf.py +++ b/PQAnalysis/analysis/rdf/rdf.py @@ -237,6 +237,11 @@ def __init__( # pylint: disable=too-many-branches RDFError If n_bins, delta_r and r_max are all specified. This would lead to ambiguous results. + RDFError + If delta_r is specified but not greater than zero. + RDFError + If the reference or target selection does not + select any atoms. Notes ----- @@ -352,10 +357,22 @@ def __init__( # pylint: disable=too-many-branches self.topology, self.use_full_atom_info ) + if len(self.reference_indices) == 0: + self.logger.error( + "The reference selection does not select any atoms.", + exception=RDFError + ) + self.target_indices = self.target_selection.select( self.topology, self.use_full_atom_info ) + if len(self.target_indices) == 0: + self.logger.error( + "The target selection does not select any atoms.", + exception=RDFError + ) + def _use_raw_fast_path(self, traj: Trajectory | TrajectoryReader) -> bool: """ Whether the raw-frame fast path is used for the given @@ -436,6 +453,15 @@ def _setup_analysis_bins( r_min: PositiveReal, ): """Selects legacy-compatible or general RDF bin semantics.""" + if delta_r is not None and delta_r <= 0.0: + self.logger.error( + ( + "The delta_r value of the RDF analysis has to be " + f"greater than zero - it actually is {delta_r}!" + ), + exception=RDFError + ) + self._legacy_rdf = self._use_legacy_rdf( n_bins=n_bins, delta_r=delta_r, @@ -824,11 +850,28 @@ def _initialize_run(self): This method is called by the run method of the RDF class. It initializes the RDF analysis for running by calculating - the average volume of the trajectory, the reference - density of the RDF analysis and the target index + the average volume of the trajectory, the reference + density of the RDF analysis and the target index combinations of the RDF analysis. + + Raises + ------ + RDFError + If the trajectory is in vacuum, as the normalization + of g(r) requires a finite cell volume. """ + if check_trajectory_vacuum(self._setup_cells): + self.logger.error( + ( + "The provided trajectory is in vacuum, so the " + "normalization of the RDF analysis requires a " + "finite cell volume. Please provide a trajectory " + "with box information." + ), + exception=RDFError + ) + self._average_volume = self._calculate_average_volume() _ref_indices_len = len(self.reference_indices) diff --git a/docs/source/userGuide/userGuide.rst b/docs/source/userGuide/userGuide.rst index e3ed1575..9b5a8b7e 100644 --- a/docs/source/userGuide/userGuide.rst +++ b/docs/source/userGuide/userGuide.rst @@ -47,8 +47,9 @@ For a file-backed periodic orthorhombic trajectory, this minimal form with RDF path. Coordinates are parsed directly as float64, while ``delta_r`` is represented as float32 as it was by the legacy C input reader. Histogram binning and all five output columns preserve the legacy arithmetic order. -Explicit ``r_max`` or ``n_bins`` values, triclinic cells, vacuum trajectories -and intra-molecular exclusion use the general PQAnalysis RDF definition. +Explicit ``r_max`` or ``n_bins`` values, triclinic cells and intra-molecular +exclusion use the general PQAnalysis RDF definition. Vacuum trajectories are +rejected, as the normalization of g(r) requires a finite cell volume. Restart files and moldescriptor files are only needed when the calculation requires molecular topology information. For example, diff --git a/tests/analysis/rdf/test_rdf.py b/tests/analysis/rdf/test_rdf.py index 09706041..ebace067 100644 --- a/tests/analysis/rdf/test_rdf.py +++ b/tests/analysis/rdf/test_rdf.py @@ -490,8 +490,16 @@ def test__init__(self, caplog): r_min=3.0, ) - system1 = AtomicSystem(cell=Cell(10, 10, 10, 90, 90, 90)) - system2 = AtomicSystem(cell=Cell(16, 13, 12, 90, 90, 90)) + system1 = AtomicSystem( + atoms=[Atom("h")], + pos=np.array([[0, 0, 0]]), + cell=Cell(10, 10, 10, 90, 90, 90) + ) + system2 = AtomicSystem( + atoms=[Atom("h")], + pos=np.array([[0, 0, 0]]), + cell=Cell(16, 13, 12, 90, 90, 90) + ) traj = Trajectory([system1, system2]) @@ -588,8 +596,12 @@ def test__init__(self, caplog): assert np.isclose(rdf.r_max, 5.0) - system1 = AtomicSystem(cell=Cell()) - system2 = AtomicSystem(cell=Cell()) + system1 = AtomicSystem( + atoms=[Atom("h")], pos=np.array([[0, 0, 0]]), cell=Cell() + ) + system2 = AtomicSystem( + atoms=[Atom("h")], pos=np.array([[0, 0, 0]]), cell=Cell() + ) traj = Trajectory([system1, system2]) @@ -663,6 +675,135 @@ def test__init__(self, caplog): assert rdf.n_bins == 5 assert np.isclose(rdf.delta_r, 1.0) + def test__init__empty_selection(self, caplog): + system = AtomicSystem( + atoms=[Atom("O"), Atom("H")], + pos=np.array([[0, 0, 0], [1, 0, 0]]), + cell=Cell(10, 10, 10, 90, 90, 90) + ) + traj = Trajectory([system]) + + assert_logging_with_exception( + caplog=caplog, + logging_name=RDF.__qualname__, + logging_level="ERROR", + message_to_test=( + "The reference selection does not select any atoms." + ), + exception=RDFError, + function=RDF, + traj=traj, + reference_species=["Na"], + target_species=["H"], + delta_r=0.5, + r_max=4.0, + ) + + assert_logging_with_exception( + caplog=caplog, + logging_name=RDF.__qualname__, + logging_level="ERROR", + message_to_test=( + "The target selection does not select any atoms." + ), + exception=RDFError, + function=RDF, + traj=traj, + reference_species=["O"], + target_species=["Na"], + delta_r=0.5, + r_max=4.0, + ) + + @pytest.mark.parametrize("example_dir", ["rdf"], indirect=False) + def test__init__empty_selection_legacy_path( + self, caplog, test_with_data_dir + ): + assert_logging_with_exception( + caplog=caplog, + logging_name=RDF.__qualname__, + logging_level="ERROR", + message_to_test=( + "The reference selection does not select any atoms." + ), + exception=RDFError, + function=RDF, + traj=TrajectoryReader("traj.xyz"), + reference_species=["Na"], + target_species=["X"], + delta_r=0.1, + ) + + def test__init__delta_r_zero(self, caplog): + system = AtomicSystem( + atoms=[Atom("h")], + pos=np.array([[0, 0, 0]]), + cell=Cell(10, 10, 10, 90, 90, 90) + ) + traj = Trajectory([system]) + + assert_logging_with_exception( + caplog=caplog, + logging_name=RDF.__qualname__, + logging_level="ERROR", + message_to_test=( + "The delta_r value of the RDF analysis has to be " + "greater than zero - it actually is 0.0!" + ), + exception=RDFError, + function=RDF, + traj=traj, + reference_species=["h"], + target_species=["h"], + delta_r=0.0, + ) + + assert_logging_with_exception( + caplog=caplog, + logging_name=RDF.__qualname__, + logging_level="ERROR", + message_to_test=( + "The delta_r value of the RDF analysis has to be " + "greater than zero - it actually is 0.0!" + ), + exception=RDFError, + function=RDF, + traj=traj, + reference_species=["h"], + target_species=["h"], + n_bins=5, + delta_r=0.0, + ) + + def test_run_vacuum_trajectory(self, caplog): + system1 = AtomicSystem( + atoms=[Atom("h"), Atom("h")], + pos=np.array([[0, 0, 0], [1, 0, 0]]), + cell=Cell() + ) + system2 = AtomicSystem( + atoms=[Atom("h"), Atom("h")], + pos=np.array([[0, 0, 0], [1, 0, 0]]), + cell=Cell() + ) + traj = Trajectory([system1, system2]) + + rdf = RDF(traj, ["h"], ["h"], delta_r=1.0, r_max=5.0) + + assert_logging_with_exception( + caplog=caplog, + logging_name=RDF.__qualname__, + logging_level="ERROR", + message_to_test=( + "The provided trajectory is in vacuum, so the " + "normalization of the RDF analysis requires a " + "finite cell volume. Please provide a trajectory " + "with box information." + ), + exception=RDFError, + function=rdf.run, + ) + @pytest.mark.parametrize("example_dir", ["rdf"], indirect=False) def test__init__uses_first_frame_topology_without_reader_topology( self, test_with_data_dir