diff --git a/export/tests/unit/test-1body.cc b/export/tests/unit/test-1body.cc
index 32c6c67be..caf76672d 100644
--- a/export/tests/unit/test-1body.cc
+++ b/export/tests/unit/test-1body.cc
@@ -695,6 +695,122 @@ TEST_CASE_METHOD(libint2::unit::DefaultFixture,
 #endif
 }
 
+TEST_CASE("q_gau input validation", "[engine][1-body][validation]") {
+  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("negative infinity exponent throws") {
+    REQUIRE_THROWS_AS(validate_gaussian_potential_primitive(
+                          {-std::numeric_limits<double>::infinity(), 1.0}),
+                      std::invalid_argument);
+  }
+
+  SECTION("non-finite coefficient throws") {
+    REQUIRE_THROWS_AS(validate_gaussian_potential_primitive(
+                          {1.0, std::numeric_limits<double>::quiet_NaN()}),
+                      std::invalid_argument);
+    REQUIRE_THROWS_AS(validate_gaussian_potential_primitive(
+                          {1.0, std::numeric_limits<double>::infinity()}),
+                      std::invalid_argument);
+    REQUIRE_THROWS_AS(validate_gaussian_potential_primitive(
+                          {1.0, -std::numeric_limits<double>::infinity()}),
+                      std::invalid_argument);
+  }
+
+  SECTION("infinite_exponent equals IEEE 754 positive infinity") {
+    REQUIRE(infinite_exponent == std::numeric_limits<double>::infinity());
+    REQUIRE(std::isinf(infinite_exponent));
+    REQUIRE(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(
+                          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(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]") {
diff --git a/include/libint2/basis.h.in b/include/libint2/basis.h.in
index b35f72e9d..7019fbe1c 100644
--- a/include/libint2/basis.h.in
+++ b/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);
 
+    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 +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::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 +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::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 +928,21 @@ 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 (pure erfc) and sigma (pure -erf).
   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)}});
diff --git a/include/libint2/boys.h b/include/libint2/boys.h
index ac2d684f8..c54a72f80 100644
--- a/include/libint2/boys.h
+++ b/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,6 +2199,13 @@ 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) {
+      assert(!isnan(prim.exponent) &&
+             "q_gau_gm_eval: primitive exponent is NaN");
+      assert(prim.exponent >= 0.0 &&
+             "q_gau_gm_eval: primitive exponent is negative");
+      assert(isfinite(prim.coefficient) &&
+             "q_gau_gm_eval: primitive coefficient is not finite");
+
       if (isinf(prim.exponent)) {
         // α = ∞ => (∞/(∞+ρ)) = 1, contributes c_i * F_m(T)
         fm_eval_->eval(&base_type::Fm_[0], T, mmax);
diff --git a/include/libint2/shell.h b/include/libint2/shell.h
index aac30a1a7..d555bd4f0 100644
--- a/include/libint2/shell.h
+++ b/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,25 @@ 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 (including
+///         -infinity), or if coefficient is non-finite (NaN or +/-infinity).
+///         Only +infinity (i.e. libint2::infinite_exponent) is accepted as a
+///         non-finite exponent.
+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 (prim.exponent < 0.0)
+    throw std::invalid_argument(
+        "GaussianPotentialPrimitive: exponent is negative");
+  if (!isfinite(prim.coefficient))
+    throw std::invalid_argument(
+        "GaussianPotentialPrimitive: coefficient is not finite");
+}
+
 /// 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:
@@ -1357,19 +1384,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.
+/// When constructing primitives by hand, call
+/// validate_gaussian_potential_primitive() on each to catch invalid inputs:
 /// @code
 ///   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{{libint2::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{{libint2::infinite_exponent, 1.0}});
 ///   centers[2] = pt;
 ///   // Ghost atom — no potential
 ///   centers[3] = nullptr;