Failproof infinity handling

export/tests/unit/test-1body.cc

@@ -695,6 +695,93 @@ TEST_CASE_METHOD(libint2::unit::DefaultFixture,
 #endif
 }
 
+TEST_CASE("q_gau input validation", "[engine][1-body][validation]") {
+  using libint2::GaussianPotentialData;
+  using libint2::GaussianPotentialPrimitive;
+  using libint2::infinite_exponent;
+  using libint2::validate_gaussian_potential_primitive;
+
+  SECTION("valid primitives do not throw") {
+    REQUIRE_NOTHROW(validate_gaussian_potential_primitive({0.0, 1.0}));
+    REQUIRE_NOTHROW(validate_gaussian_potential_primitive({1.5, -0.3}));
+    REQUIRE_NOTHROW(
+        validate_gaussian_potential_primitive({infinite_exponent, 1.0}));
+    REQUIRE_NOTHROW(
+        validate_gaussian_potential_primitive({infinite_exponent, 0.0}));
+  }
+
+  SECTION("NaN exponent throws") {
+    REQUIRE_THROWS_AS(validate_gaussian_potential_primitive(
+                          {std::numeric_limits<double>::quiet_NaN(), 1.0}),
+                      std::invalid_argument);
+  }
+
+  SECTION("negative exponent throws") {
+    REQUIRE_THROWS_AS(validate_gaussian_potential_primitive({-1.0, 1.0}),
+                      std::invalid_argument);
+  }
+
+  SECTION("NaN coefficient throws") {
+    REQUIRE_THROWS_AS(validate_gaussian_potential_primitive(
+                          {1.0, std::numeric_limits<double>::quiet_NaN()}),
+                      std::invalid_argument);
+  }
+
+  SECTION("infinite_exponent equals IEEE 754 infinity") {
+    REQUIRE(infinite_exponent == std::numeric_limits<double>::infinity());
+    REQUIRE(std::isinf(infinite_exponent));
+  }
+}
+
+TEST_CASE("make_q_gau_data factory validation",
+          "[engine][1-body][validation]") {
+  using libint2::Atom;
+
+  std::vector<Atom> atoms = {Atom{1, 0.0, 0.0, 0.0}};
+
+  SECTION("make_q_gau_data_erf rejects NaN omega") {
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erf(
+                          std::numeric_limits<double>::quiet_NaN(), atoms),
+                      std::invalid_argument);
+  }
+
+  SECTION("make_q_gau_data_erf rejects non-positive omega") {
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erf(0.0, atoms),
+                      std::invalid_argument);
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erf(-1.0, atoms),
+                      std::invalid_argument);
+  }
+
+  SECTION("make_q_gau_data_erfc rejects invalid omega") {
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erfc(
+                          std::numeric_limits<double>::quiet_NaN(), atoms),
+                      std::invalid_argument);
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erfc(0.0, atoms),
+                      std::invalid_argument);
+  }
+
+  SECTION("make_q_gau_data_erfx rejects invalid parameters") {
+    REQUIRE_THROWS_AS(
+        libint2::make_q_gau_data_erfx(std::numeric_limits<double>::quiet_NaN(),
+                                      0.3, 0.7, atoms),
+        std::invalid_argument);
+    REQUIRE_THROWS_AS(
+        libint2::make_q_gau_data_erfx(
+            0.5, std::numeric_limits<double>::quiet_NaN(), 0.7, atoms),
+        std::invalid_argument);
+    REQUIRE_THROWS_AS(
+        libint2::make_q_gau_data_erfx(
+            0.5, 0.3, std::numeric_limits<double>::quiet_NaN(), atoms),
+        std::invalid_argument);
+  }
+
+  SECTION("valid factory calls succeed") {
+    REQUIRE_NOTHROW(libint2::make_q_gau_data_erf(0.5, atoms));
+    REQUIRE_NOTHROW(libint2::make_q_gau_data_erfc(0.5, atoms));
+    REQUIRE_NOTHROW(libint2::make_q_gau_data_erfx(0.5, 0.3, 0.7, atoms));
+  }
+}
+
 // verify that python/tests/test_libint2.py:test_integrals is correct
 TEST_CASE_METHOD(libint2::unit::DefaultFixture, "python correctness",
                  "[engine][1-body]") {

include/libint2/basis.h.in

@@ -801,7 +801,7 @@ namespace libint2 {
   /// @param[in] atoms the atoms
   inline GaussianPotentialCentersData make_q_gau_data(
       NuclearModel model, const std::vector<Atom>& atoms) {
-    constexpr double inf = std::numeric_limits<double>::infinity();
+    constexpr double inf = libint2::infinite_exponent;
     std::map<int, std::shared_ptr<const GaussianPotentialData>> cache;
     GaussianPotentialCentersData result;
     result.reserve(atoms.size());
@@ -837,7 +837,7 @@ namespace libint2 {
   inline GaussianPotentialCentersData make_q_gau_data(
       NuclearModel model, const std::vector<Atom>& atoms,
       const std::string& sap_basis_name) {
-    constexpr double inf = std::numeric_limits<double>::infinity();
+    constexpr double inf = libint2::infinite_exponent;
     std::string basis_lib_path = basis_data_path();
     std::string canonical_name = sap_basis_name;
     std::transform(sap_basis_name.begin(), sap_basis_name.end(),
@@ -849,6 +849,17 @@ namespace libint2 {
     auto file = basis_lib_path + PATH_SEPARATOR + canonical_name + ".g94";
     auto sap_by_Z = read_sap_basis_library(file);
 
+#ifndef NDEBUG
+    for (size_t z = 0; z < sap_by_Z.size(); ++z) {
+      for (const auto& p : sap_by_Z[z]) {
+        assert(!std::isnan(p.exponent) && p.exponent > 0.0 &&
+               "SAP basis data contains invalid exponent");
+        assert(!std::isnan(p.coefficient) &&
+               "SAP basis data contains NaN coefficient");
+      }
+    }
+#endif
+
     std::map<int, std::shared_ptr<const GaussianPotentialData>> cache;
     GaussianPotentialCentersData result;
     result.reserve(atoms.size());
@@ -888,6 +899,9 @@ namespace libint2 {
   /// Build GaussianPotentialCentersData for erf(ω)/r model.
   inline GaussianPotentialCentersData make_q_gau_data_erf(
       double omega, const std::vector<Atom>& atoms) {
+    if (std::isnan(omega) || omega <= 0.0)
+      throw std::invalid_argument(
+          "make_q_gau_data_erf: omega must be positive and finite");
     auto data = std::make_shared<const GaussianPotentialData>(
         GaussianPotentialData{{omega * omega, 1.0}});
     GaussianPotentialCentersData result;
@@ -900,7 +914,10 @@ namespace libint2 {
   /// Build GaussianPotentialCentersData for erfc(ω)/r model.
   inline GaussianPotentialCentersData make_q_gau_data_erfc(
       double omega, const std::vector<Atom>& atoms) {
-    constexpr double inf = std::numeric_limits<double>::infinity();
+    if (std::isnan(omega) || omega <= 0.0)
+      throw std::invalid_argument(
+          "make_q_gau_data_erfc: omega must be positive and finite");
+    constexpr double inf = libint2::infinite_exponent;
     auto data = std::make_shared<const GaussianPotentialData>(
         GaussianPotentialData{{inf, 1.0}, {omega * omega, -1.0}});
     GaussianPotentialCentersData result;
@@ -915,7 +932,16 @@ namespace libint2 {
   inline GaussianPotentialCentersData make_q_gau_data_erfx(
       double omega, double lambda, double sigma,
       const std::vector<Atom>& atoms) {
-    constexpr double inf = std::numeric_limits<double>::infinity();
+    if (std::isnan(omega) || omega <= 0.0)
+      throw std::invalid_argument(
+          "make_q_gau_data_erfx: omega must be positive and finite");
+    if (std::isnan(lambda))
+      throw std::invalid_argument(
+          "make_q_gau_data_erfx: lambda must not be NaN");
+    if (std::isnan(sigma))
+      throw std::invalid_argument(
+          "make_q_gau_data_erfx: sigma must not be NaN");
+    constexpr double inf = libint2::infinite_exponent;
     auto data = std::make_shared<const GaussianPotentialData>(
         GaussianPotentialData{{inf, sigma},
                               {omega * omega, -(sigma - lambda)}});

include/libint2/boys.h

@@ -2196,6 +2196,18 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch<q_gau_gm_eval<Real>> {
 
     std::fill(Gm, Gm + mmax + 1, Real{0});
 
+#ifndef NDEBUG
+    for (const auto& prim : primitives) {
+      using std::isnan;
+      assert(!isnan(prim.exponent) &&
+             "q_gau_gm_eval: primitive exponent is NaN");
+      assert(!isnan(prim.coefficient) &&
+             "q_gau_gm_eval: primitive coefficient is NaN");
+      assert((isinf(prim.exponent) || prim.exponent >= 0.0) &&
+             "q_gau_gm_eval: primitive exponent is negative");
+    }
+#endif
+
     for (const auto& prim : primitives) {
       if (isinf(prim.exponent)) {
         // α = ∞ => (∞/(∞+ρ)) = 1, contributes c_i * F_m(T)

include/libint2/shell.h

@@ -34,8 +34,10 @@
 #include <cassert>
 #include <cmath>
 #include <iostream>
+#include <limits>
 #include <map>
 #include <memory>
+#include <stdexcept>
 #include <type_traits>
 #include <vector>
 
@@ -1332,12 +1334,18 @@ struct ShellPair {
   }
 };
 
+/// Sentinel value for GaussianPotentialPrimitive::exponent indicating a point
+/// charge (bare Coulomb 1/r potential).  Prefer this over raw INFINITY or
+/// std::numeric_limits<double>::infinity() for clarity.
+constexpr double infinite_exponent = std::numeric_limits<double>::infinity();
+
 /// A single Gaussian primitive contributing to a nuclear potential.
 /// Each primitive represents a term c * (α/(α+ρ))^(m+1/2) * F_m(T·α/(α+ρ))
-/// in the core integral expansion. Use exponent = infinity for point-charge
-/// contributions (recovers bare Coulomb F_m(T)).
+/// in the core integral expansion. Use exponent = libint2::infinite_exponent
+/// for point-charge contributions (recovers bare Coulomb F_m(T)).
 struct GaussianPotentialPrimitive {
-  double exponent;     ///< Gaussian exponent α (infinity for point charge)
+  double exponent;  ///< Gaussian exponent α; use libint2::infinite_exponent for
+                    ///< point charge
   double coefficient;  ///< coefficient c
 };
 
@@ -1346,6 +1354,23 @@ struct GaussianPotentialPrimitive {
 /// corrections, erf/erfc attenuations, or any combination thereof.
 using GaussianPotentialData = std::vector<GaussianPotentialPrimitive>;
 
+/// Validates a GaussianPotentialPrimitive.
+/// @throws std::invalid_argument if exponent is NaN or negative-finite,
+///         or if coefficient is NaN
+inline void validate_gaussian_potential_primitive(
+    const GaussianPotentialPrimitive& prim) {
+  using std::isfinite;
+  using std::isnan;
+  if (isnan(prim.exponent))
+    throw std::invalid_argument("GaussianPotentialPrimitive: exponent is NaN");
+  if (isfinite(prim.exponent) && prim.exponent < 0.0)
+    throw std::invalid_argument(
+        "GaussianPotentialPrimitive: exponent is negative");
+  if (isnan(prim.coefficient))
+    throw std::invalid_argument(
+        "GaussianPotentialPrimitive: coefficient is NaN");
+}
+
 /// Gaussian potential data per center, parallel to the point-charges vector
 /// passed to Engine::set_params(). Each element is a shared_ptr to the
 /// potential primitives for that center:
@@ -1362,14 +1387,15 @@ using GaussianPotentialData = std::vector<GaussianPotentialPrimitive>;
 ///   GaussianPotentialCentersData centers(point_charges.size());
 ///   // QM atom — point nuclear + SAP correction
 ///   centers[0] = std::make_shared<const GaussianPotentialData>(
-///       GaussianPotentialData{{INFINITY, 1.0}, {alpha1, c1}, {alpha2, c2}});
+///       GaussianPotentialData{{infinite_exponent, 1.0}, {alpha1, c1}, {alpha2,
+///       c2}});
 ///   // QM atom — finite (Gaussian) nuclear model
 ///   auto xi = chemistry::gaussian_nuclear_exponent(Z);
 ///   centers[1] = std::make_shared<const GaussianPotentialData>(
 ///       GaussianPotentialData{{xi, 1.0}});
 ///   // MM point charge — bare Coulomb (point nuclear)
 ///   static auto pt = std::make_shared<const GaussianPotentialData>(
-///       GaussianPotentialData{{INFINITY, 1.0}});
+///       GaussianPotentialData{{infinite_exponent, 1.0}});
 ///   centers[2] = pt;
 ///   // Ghost atom — no potential
 ///   centers[3] = nullptr;