Failproof infinity handling

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

@@ -695,6 +695,86 @@ TEST_CASE_METHOD(libint2::unit::DefaultFixture,
 #endif
 }
 
+TEST_CASE("infinite_exponent sentinel", "[engine][1-body][validation]") {
+  REQUIRE(libint2::infinite_exponent ==
+          std::numeric_limits<double>::infinity());
+  REQUIRE(std::isinf(libint2::infinite_exponent));
+  REQUIRE(libint2::infinite_exponent > 0.0);
+}
+
+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 or non-finite 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);
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erf(
+                          std::numeric_limits<double>::infinity(), 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);
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erfc(-1.0, atoms),
+                      std::invalid_argument);
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erfc(
+                          std::numeric_limits<double>::infinity(), 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.0, 0.3, 0.7, atoms),
+                      std::invalid_argument);
+    REQUIRE_THROWS_AS(libint2::make_q_gau_data_erfx(-1.0, 0.3, 0.7, atoms),
+                      std::invalid_argument);
+    REQUIRE_THROWS_AS(
+        libint2::make_q_gau_data_erfx(std::numeric_limits<double>::infinity(),
+                                      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, std::numeric_limits<double>::infinity(), 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);
+    REQUIRE_THROWS_AS(
+        libint2::make_q_gau_data_erfx(
+            0.5, 0.3, std::numeric_limits<double>::infinity(), 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,20 @@ namespace libint2 {
     auto file = basis_lib_path + PATH_SEPARATOR + canonical_name + ".g94";
     auto sap_by_Z = read_sap_basis_library(file);
 
+    // Validate unconditionally: SAP data comes from disk (system boundary).
+    // A corrupt .g94 file must not silently propagate invalid exponents or
+    // coefficients into integral evaluation.
+    for (size_t z = 0; z < sap_by_Z.size(); ++z) {
+      for (const auto& p : sap_by_Z[z]) {
+        if (!std::isfinite(p.exponent) || p.exponent <= 0.0)
+          throw std::invalid_argument(
+              "make_q_gau_data: SAP basis contains invalid exponent");
+        if (!std::isfinite(p.coefficient))
+          throw std::invalid_argument(
+              "make_q_gau_data: SAP basis contains non-finite coefficient");
+      }
+    }
+
     std::map<int, std::shared_ptr<const GaussianPotentialData>> cache;
     GaussianPotentialCentersData result;
     result.reserve(atoms.size());
@@ -888,6 +902,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::isfinite(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 +917,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::isfinite(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;
@@ -911,11 +931,22 @@ namespace libint2 {
   }
 
   /// Build GaussianPotentialCentersData for erfx(ω,λ,σ)/r model.
-  /// Computes (λ·erf(ωr) + σ·erfc(ωr))/r = σ/r - (σ-λ)·erf(ωr)/r
+  /// Computes (λ·erf(ωr) + σ·erfc(ωr))/r = σ/r - (σ-λ)·erf(ωr)/r.
+  /// Zero is valid for lambda and sigma: λ = 0 gives σ·erfc(ωr)/r,
+  /// σ = 0 gives λ·erf(ωr)/r.
   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::isfinite(omega) || omega <= 0.0)
+      throw std::invalid_argument(
+          "make_q_gau_data_erfx: omega must be positive and finite");
+    if (!std::isfinite(lambda))
+      throw std::invalid_argument(
+          "make_q_gau_data_erfx: lambda must be finite");
+    if (!std::isfinite(sigma))
+      throw std::invalid_argument(
+          "make_q_gau_data_erfx: sigma must be finite");
+    constexpr double inf = libint2::infinite_exponent;
     auto data = std::make_shared<const GaussianPotentialData>(
         GaussianPotentialData{{inf, sigma},
                               {omega * omega, -(sigma - lambda)}});

include/libint2/boys.h

@@ -2186,7 +2186,9 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch<q_gau_gm_eval<Real>> {
   template <typename PrimitivesContainer>
   void operator()(Real* Gm, Real rho, Real T, int mmax,
                   const PrimitivesContainer& primitives) {
+    using std::isfinite;
     using std::isinf;
+    using std::isnan;
     using std::sqrt;
 
     if (primitives.empty()) {
@@ -2197,7 +2199,17 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch<q_gau_gm_eval<Real>> {
     std::fill(Gm, Gm + mmax + 1, Real{0});
 
     for (const auto& prim : primitives) {
-      if (isinf(prim.exponent)) {
+      assert(!isnan(prim.exponent) &&
+             "q_gau_gm_eval: primitive exponent is NaN");
+      assert(prim.exponent > 0.0 &&
+             "q_gau_gm_eval: primitive exponent is non-positive");
+      assert(isfinite(prim.coefficient) &&
+             "q_gau_gm_eval: primitive coefficient is not finite");
+
+      // Only +inf is the point-charge sentinel. -inf (invalid input that
+      // escapes validation in release) falls into the Gaussian branch and
+      // produces NaN from -inf/(-inf+rho) — visibly wrong rather than silently.
+      if (isinf(prim.exponent) && prim.exponent > 0.0) {
         // α = ∞ => (∞/(∞+ρ)) = 1, contributes c_i * F_m(T)
         fm_eval_->eval(&base_type::Fm_[0], T, mmax);
         for (auto m = 0; m <= mmax; ++m)

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
 };
 
@@ -1357,19 +1365,22 @@ using GaussianPotentialData = std::vector<GaussianPotentialPrimitive>;
 ///
 /// The convenience functions make_q_gau_data() build this from a nuclear model
 /// specification and a list of atoms. For full per-center control (e.g.,
-/// different models for QM vs MM atoms), construct the vector manually:
+/// different models for QM vs MM atoms), construct the vector manually.
+/// Each primitive must satisfy: exponent > 0 or exponent ==
+/// libint2::infinite_exponent; coefficient must be finite (not NaN or +/-inf).
 /// @code
+///   using libint2::infinite_exponent;
 ///   GaussianPotentialCentersData centers(point_charges.size());
-///   // QM atom — point nuclear + SAP correction
+///   // QM atom — point nuclear + one SAP primitive
 ///   centers[0] = std::make_shared<const GaussianPotentialData>(
-///       GaussianPotentialData{{INFINITY, 1.0}, {alpha1, c1}, {alpha2, c2}});
+///       GaussianPotentialData{{infinite_exponent, 1.0}, {alpha1, c1}});
 ///   // 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;