From 0424882b98e3217a363f0bb7bfa37063a9893ae8 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 25 Mar 2026 11:44:09 -0400 Subject: [PATCH 1/9] makes infinity handling failproof for Operator::q_gau --- export/tests/unit/test-1body.cc | 87 +++++++++++++++++++++++++++++++++ include/libint2/basis.h.in | 34 +++++++++++-- include/libint2/boys.h | 12 +++++ include/libint2/shell.h | 36 ++++++++++++-- 4 files changed, 160 insertions(+), 9 deletions(-) diff --git a/export/tests/unit/test-1body.cc b/export/tests/unit/test-1body.cc index 32c6c67be..0b5e489b3 100644 --- a/export/tests/unit/test-1body.cc +++ b/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::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::quiet_NaN()}), + std::invalid_argument); + } + + SECTION("infinite_exponent equals IEEE 754 infinity") { + REQUIRE(infinite_exponent == std::numeric_limits::infinity()); + REQUIRE(std::isinf(infinite_exponent)); + } +} + +TEST_CASE("make_q_gau_data factory validation", + "[engine][1-body][validation]") { + using libint2::Atom; + + std::vector 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::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::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::quiet_NaN(), + 0.3, 0.7, atoms), + std::invalid_argument); + REQUIRE_THROWS_AS( + libint2::make_q_gau_data_erfx( + 0.5, std::numeric_limits::quiet_NaN(), 0.7, atoms), + std::invalid_argument); + REQUIRE_THROWS_AS( + libint2::make_q_gau_data_erfx( + 0.5, 0.3, std::numeric_limits::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]") { diff --git a/include/libint2/basis.h.in b/include/libint2/basis.h.in index b35f72e9d..05e8c9b35 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& atoms) { - constexpr double inf = std::numeric_limits::infinity(); + constexpr double inf = libint2::infinite_exponent; std::map> cache; GaussianPotentialCentersData result; result.reserve(atoms.size()); @@ -837,7 +837,7 @@ namespace libint2 { inline GaussianPotentialCentersData make_q_gau_data( NuclearModel model, const std::vector& atoms, const std::string& sap_basis_name) { - constexpr double inf = std::numeric_limits::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> 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& 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( 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& atoms) { - constexpr double inf = std::numeric_limits::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( 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& atoms) { - constexpr double inf = std::numeric_limits::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( GaussianPotentialData{{inf, sigma}, {omega * omega, -(sigma - lambda)}}); diff --git a/include/libint2/boys.h b/include/libint2/boys.h index ac2d684f8..3b2cc0e32 100644 --- a/include/libint2/boys.h +++ b/include/libint2/boys.h @@ -2196,6 +2196,18 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch> { 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) diff --git a/include/libint2/shell.h b/include/libint2/shell.h index aac30a1a7..fc40421d4 100644 --- a/include/libint2/shell.h +++ b/include/libint2/shell.h @@ -34,8 +34,10 @@ #include #include #include +#include #include #include +#include #include #include @@ -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::infinity() for clarity. +constexpr double infinite_exponent = std::numeric_limits::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; +/// 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; /// GaussianPotentialCentersData centers(point_charges.size()); /// // QM atom — point nuclear + SAP correction /// centers[0] = std::make_shared( -/// 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( /// GaussianPotentialData{{xi, 1.0}}); /// // MM point charge — bare Coulomb (point nuclear) /// static auto pt = std::make_shared( -/// GaussianPotentialData{{INFINITY, 1.0}}); +/// GaussianPotentialData{{infinite_exponent, 1.0}}); /// centers[2] = pt; /// // Ghost atom — no potential /// centers[3] = nullptr; From 4d9a35e2509a8d1966438a178c526c25ac2235c3 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 25 Mar 2026 14:28:10 -0400 Subject: [PATCH 2/9] reject -inf exponents, non-finite coefficients/params --- export/tests/unit/test-1body.cc | 34 ++++++++++++++++++++++++++++++--- include/libint2/basis.h.in | 20 +++++++++---------- include/libint2/boys.h | 7 ++++--- include/libint2/shell.h | 12 +++++++----- 4 files changed, 52 insertions(+), 21 deletions(-) diff --git a/export/tests/unit/test-1body.cc b/export/tests/unit/test-1body.cc index 0b5e489b3..fe039a979 100644 --- a/export/tests/unit/test-1body.cc +++ b/export/tests/unit/test-1body.cc @@ -721,15 +721,28 @@ TEST_CASE("q_gau input validation", "[engine][1-body][validation]") { std::invalid_argument); } - SECTION("NaN coefficient throws") { + SECTION("negative infinity exponent throws") { + REQUIRE_THROWS_AS(validate_gaussian_potential_primitive( + {-std::numeric_limits::infinity(), 1.0}), + std::invalid_argument); + } + + SECTION("non-finite coefficient throws") { REQUIRE_THROWS_AS(validate_gaussian_potential_primitive( {1.0, std::numeric_limits::quiet_NaN()}), std::invalid_argument); + REQUIRE_THROWS_AS(validate_gaussian_potential_primitive( + {1.0, std::numeric_limits::infinity()}), + std::invalid_argument); + REQUIRE_THROWS_AS(validate_gaussian_potential_primitive( + {1.0, -std::numeric_limits::infinity()}), + std::invalid_argument); } - SECTION("infinite_exponent equals IEEE 754 infinity") { + SECTION("infinite_exponent equals IEEE 754 positive infinity") { REQUIRE(infinite_exponent == std::numeric_limits::infinity()); REQUIRE(std::isinf(infinite_exponent)); + REQUIRE(infinite_exponent > 0.0); } } @@ -745,11 +758,14 @@ TEST_CASE("make_q_gau_data factory validation", std::invalid_argument); } - SECTION("make_q_gau_data_erf rejects non-positive omega") { + 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::infinity(), atoms), + std::invalid_argument); } SECTION("make_q_gau_data_erfc rejects invalid omega") { @@ -765,14 +781,26 @@ TEST_CASE("make_q_gau_data factory validation", libint2::make_q_gau_data_erfx(std::numeric_limits::quiet_NaN(), 0.3, 0.7, atoms), std::invalid_argument); + REQUIRE_THROWS_AS( + libint2::make_q_gau_data_erfx(std::numeric_limits::infinity(), + 0.3, 0.7, atoms), + std::invalid_argument); REQUIRE_THROWS_AS( libint2::make_q_gau_data_erfx( 0.5, std::numeric_limits::quiet_NaN(), 0.7, atoms), std::invalid_argument); + REQUIRE_THROWS_AS( + libint2::make_q_gau_data_erfx( + 0.5, std::numeric_limits::infinity(), 0.7, atoms), + std::invalid_argument); REQUIRE_THROWS_AS( libint2::make_q_gau_data_erfx( 0.5, 0.3, std::numeric_limits::quiet_NaN(), atoms), std::invalid_argument); + REQUIRE_THROWS_AS( + libint2::make_q_gau_data_erfx( + 0.5, 0.3, std::numeric_limits::infinity(), atoms), + std::invalid_argument); } SECTION("valid factory calls succeed") { diff --git a/include/libint2/basis.h.in b/include/libint2/basis.h.in index 05e8c9b35..955faa496 100644 --- a/include/libint2/basis.h.in +++ b/include/libint2/basis.h.in @@ -852,10 +852,10 @@ namespace libint2 { #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 && + assert(std::isfinite(p.exponent) && p.exponent > 0.0 && "SAP basis data contains invalid exponent"); - assert(!std::isnan(p.coefficient) && - "SAP basis data contains NaN coefficient"); + assert(std::isfinite(p.coefficient) && + "SAP basis data contains non-finite coefficient"); } } #endif @@ -899,7 +899,7 @@ namespace libint2 { /// Build GaussianPotentialCentersData for erf(ω)/r model. inline GaussianPotentialCentersData make_q_gau_data_erf( double omega, const std::vector& atoms) { - if (std::isnan(omega) || omega <= 0.0) + 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( @@ -914,7 +914,7 @@ namespace libint2 { /// Build GaussianPotentialCentersData for erfc(ω)/r model. inline GaussianPotentialCentersData make_q_gau_data_erfc( double omega, const std::vector& atoms) { - if (std::isnan(omega) || omega <= 0.0) + 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; @@ -932,15 +932,15 @@ namespace libint2 { inline GaussianPotentialCentersData make_q_gau_data_erfx( double omega, double lambda, double sigma, const std::vector& atoms) { - if (std::isnan(omega) || omega <= 0.0) + if (!std::isfinite(omega) || omega <= 0.0) throw std::invalid_argument( "make_q_gau_data_erfx: omega must be positive and finite"); - if (std::isnan(lambda)) + if (!std::isfinite(lambda)) throw std::invalid_argument( - "make_q_gau_data_erfx: lambda must not be NaN"); - if (std::isnan(sigma)) + "make_q_gau_data_erfx: lambda must be finite"); + if (!std::isfinite(sigma)) throw std::invalid_argument( - "make_q_gau_data_erfx: sigma must not be NaN"); + "make_q_gau_data_erfx: sigma must be finite"); constexpr double inf = libint2::infinite_exponent; auto data = std::make_shared( GaussianPotentialData{{inf, sigma}, diff --git a/include/libint2/boys.h b/include/libint2/boys.h index 3b2cc0e32..857f7b58c 100644 --- a/include/libint2/boys.h +++ b/include/libint2/boys.h @@ -2198,13 +2198,14 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch> { #ifndef NDEBUG for (const auto& prim : primitives) { + using std::isfinite; 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) && + 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"); } #endif diff --git a/include/libint2/shell.h b/include/libint2/shell.h index fc40421d4..22b8634b3 100644 --- a/include/libint2/shell.h +++ b/include/libint2/shell.h @@ -1355,20 +1355,22 @@ struct GaussianPotentialPrimitive { using GaussianPotentialData = std::vector; /// Validates a GaussianPotentialPrimitive. -/// @throws std::invalid_argument if exponent is NaN or negative-finite, -/// or if coefficient is NaN +/// @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 (isfinite(prim.exponent) && prim.exponent < 0.0) + if (prim.exponent < 0.0) throw std::invalid_argument( "GaussianPotentialPrimitive: exponent is negative"); - if (isnan(prim.coefficient)) + if (!isfinite(prim.coefficient)) throw std::invalid_argument( - "GaussianPotentialPrimitive: coefficient is NaN"); + "GaussianPotentialPrimitive: coefficient is not finite"); } /// Gaussian potential data per center, parallel to the point-charges vector From 11e32b4ccf6969b303664b7e53acbce1ed71475d Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 25 Mar 2026 14:52:24 -0400 Subject: [PATCH 3/9] cleanup --- export/tests/unit/test-1body.cc | 2 -- include/libint2/shell.h | 6 +++--- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/export/tests/unit/test-1body.cc b/export/tests/unit/test-1body.cc index fe039a979..2114e812f 100644 --- a/export/tests/unit/test-1body.cc +++ b/export/tests/unit/test-1body.cc @@ -696,8 +696,6 @@ TEST_CASE_METHOD(libint2::unit::DefaultFixture, } 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; diff --git a/include/libint2/shell.h b/include/libint2/shell.h index 22b8634b3..977ba8b2f 100644 --- a/include/libint2/shell.h +++ b/include/libint2/shell.h @@ -1389,15 +1389,15 @@ inline void validate_gaussian_potential_primitive( /// GaussianPotentialCentersData centers(point_charges.size()); /// // QM atom — point nuclear + SAP correction /// centers[0] = std::make_shared( -/// GaussianPotentialData{{infinite_exponent, 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( /// GaussianPotentialData{{xi, 1.0}}); /// // MM point charge — bare Coulomb (point nuclear) /// static auto pt = std::make_shared( -/// GaussianPotentialData{{infinite_exponent, 1.0}}); +/// GaussianPotentialData{{libint2::infinite_exponent, 1.0}}); /// centers[2] = pt; /// // Ghost atom — no potential /// centers[3] = nullptr; From ef560b4fbdac97f93ef76150328cd5fc3ae16396 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 25 Mar 2026 16:37:44 -0400 Subject: [PATCH 4/9] merge assertions into eval loop, unconditional SAP validation, document erfx zero params --- include/libint2/basis.h.in | 15 ++++++++------- include/libint2/boys.h | 9 ++------- 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/include/libint2/basis.h.in b/include/libint2/basis.h.in index 955faa496..7019fbe1c 100644 --- a/include/libint2/basis.h.in +++ b/include/libint2/basis.h.in @@ -849,16 +849,16 @@ 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::isfinite(p.exponent) && p.exponent > 0.0 && - "SAP basis data contains invalid exponent"); - assert(std::isfinite(p.coefficient) && - "SAP basis data contains non-finite coefficient"); + 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"); } } -#endif std::map> cache; GaussianPotentialCentersData result; @@ -928,7 +928,8 @@ 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& atoms) { diff --git a/include/libint2/boys.h b/include/libint2/boys.h index 857f7b58c..f261d49cb 100644 --- a/include/libint2/boys.h +++ b/include/libint2/boys.h @@ -2187,6 +2187,7 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch> { void operator()(Real* Gm, Real rho, Real T, int mmax, const PrimitivesContainer& primitives) { using std::isinf; + using std::isnan; using std::sqrt; if (primitives.empty()) { @@ -2196,20 +2197,14 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch> { std::fill(Gm, Gm + mmax + 1, Real{0}); -#ifndef NDEBUG for (const auto& prim : primitives) { - using std::isfinite; - using std::isnan; 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) && + assert(std::isfinite(static_cast(prim.coefficient)) && "q_gau_gm_eval: primitive coefficient is not finite"); - } -#endif - for (const auto& prim : primitives) { if (isinf(prim.exponent)) { // α = ∞ => (∞/(∞+ρ)) = 1, contributes c_i * F_m(T) fm_eval_->eval(&base_type::Fm_[0], T, mmax); From 9fac5df751a7cae267eeb3462de7a96df2fdf96f Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Thu, 26 Mar 2026 14:19:40 -0400 Subject: [PATCH 5/9] cleanup: fix static_cast in assert, add erfc +inf test, doc validate hint --- export/tests/unit/test-1body.cc | 3 +++ include/libint2/boys.h | 3 ++- include/libint2/shell.h | 4 +++- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/export/tests/unit/test-1body.cc b/export/tests/unit/test-1body.cc index 2114e812f..caf76672d 100644 --- a/export/tests/unit/test-1body.cc +++ b/export/tests/unit/test-1body.cc @@ -772,6 +772,9 @@ TEST_CASE("make_q_gau_data factory validation", 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::infinity(), atoms), + std::invalid_argument); } SECTION("make_q_gau_data_erfx rejects invalid parameters") { diff --git a/include/libint2/boys.h b/include/libint2/boys.h index f261d49cb..c54a72f80 100644 --- a/include/libint2/boys.h +++ b/include/libint2/boys.h @@ -2186,6 +2186,7 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch> { template 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; @@ -2202,7 +2203,7 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch> { "q_gau_gm_eval: primitive exponent is NaN"); assert(prim.exponent >= 0.0 && "q_gau_gm_eval: primitive exponent is negative"); - assert(std::isfinite(static_cast(prim.coefficient)) && + assert(isfinite(prim.coefficient) && "q_gau_gm_eval: primitive coefficient is not finite"); if (isinf(prim.exponent)) { diff --git a/include/libint2/shell.h b/include/libint2/shell.h index 977ba8b2f..d555bd4f0 100644 --- a/include/libint2/shell.h +++ b/include/libint2/shell.h @@ -1384,7 +1384,9 @@ inline void validate_gaussian_potential_primitive( /// /// 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 From 72ad405d0a6b8d036a244fc4c4a8d5bfe2579f1f Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Thu, 16 Apr 2026 15:24:43 -0400 Subject: [PATCH 6/9] fix erfx doc, clarify SAP validation, add missing factory tests --- export/tests/unit/test-1body.cc | 6 ++++++ include/libint2/basis.h.in | 6 +++++- include/libint2/shell.h | 8 ++++---- 3 files changed, 15 insertions(+), 5 deletions(-) diff --git a/export/tests/unit/test-1body.cc b/export/tests/unit/test-1body.cc index caf76672d..40d817e3d 100644 --- a/export/tests/unit/test-1body.cc +++ b/export/tests/unit/test-1body.cc @@ -772,6 +772,8 @@ TEST_CASE("make_q_gau_data factory validation", 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::infinity(), atoms), std::invalid_argument); @@ -782,6 +784,10 @@ TEST_CASE("make_q_gau_data factory validation", libint2::make_q_gau_data_erfx(std::numeric_limits::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::infinity(), 0.3, 0.7, atoms), diff --git a/include/libint2/basis.h.in b/include/libint2/basis.h.in index 7019fbe1c..8bc2c26ca 100644 --- a/include/libint2/basis.h.in +++ b/include/libint2/basis.h.in @@ -849,6 +849,9 @@ 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) @@ -929,7 +932,8 @@ namespace libint2 { /// Build GaussianPotentialCentersData for erfx(ω,λ,σ)/r model. /// Computes (λ·erf(ωr) + σ·erfc(ωr))/r = σ/r - (σ-λ)·erf(ωr)/r. - /// Zero is valid for lambda (pure erfc) and sigma (pure -erf). + /// 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& atoms) { diff --git a/include/libint2/shell.h b/include/libint2/shell.h index d555bd4f0..96adbeaa3 100644 --- a/include/libint2/shell.h +++ b/include/libint2/shell.h @@ -1388,18 +1388,18 @@ inline void validate_gaussian_potential_primitive( /// When constructing primitives by hand, call /// validate_gaussian_potential_primitive() on each to catch invalid inputs: /// @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( -/// GaussianPotentialData{{libint2::infinite_exponent, 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( /// GaussianPotentialData{{xi, 1.0}}); /// // MM point charge — bare Coulomb (point nuclear) /// static auto pt = std::make_shared( -/// GaussianPotentialData{{libint2::infinite_exponent, 1.0}}); +/// GaussianPotentialData{{infinite_exponent, 1.0}}); /// centers[2] = pt; /// // Ghost atom — no potential /// centers[3] = nullptr; From a5dbbfd4aaa34504e5629c1864cecf23c41a86a6 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Thu, 16 Apr 2026 15:38:43 -0400 Subject: [PATCH 7/9] dispatch only on +inf sentinel, make -inf fail loud in release --- include/libint2/boys.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/libint2/boys.h b/include/libint2/boys.h index c54a72f80..9a26e8786 100644 --- a/include/libint2/boys.h +++ b/include/libint2/boys.h @@ -2206,7 +2206,10 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch> { assert(isfinite(prim.coefficient) && "q_gau_gm_eval: primitive coefficient is not finite"); - if (isinf(prim.exponent)) { + // 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) From 612978dfa5f9bf1cfedeb23188cc2032c8245192 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Fri, 17 Apr 2026 12:33:31 -0400 Subject: [PATCH 8/9] remove unused validate helper, allow zero SAP exponent --- export/tests/unit/test-1body.cc | 52 ++++----------------------------- include/libint2/basis.h.in | 2 +- include/libint2/shell.h | 23 ++------------- 3 files changed, 8 insertions(+), 69 deletions(-) diff --git a/export/tests/unit/test-1body.cc b/export/tests/unit/test-1body.cc index 40d817e3d..8c27be49f 100644 --- a/export/tests/unit/test-1body.cc +++ b/export/tests/unit/test-1body.cc @@ -695,53 +695,11 @@ 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::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::infinity(), 1.0}), - std::invalid_argument); - } - - SECTION("non-finite coefficient throws") { - REQUIRE_THROWS_AS(validate_gaussian_potential_primitive( - {1.0, std::numeric_limits::quiet_NaN()}), - std::invalid_argument); - REQUIRE_THROWS_AS(validate_gaussian_potential_primitive( - {1.0, std::numeric_limits::infinity()}), - std::invalid_argument); - REQUIRE_THROWS_AS(validate_gaussian_potential_primitive( - {1.0, -std::numeric_limits::infinity()}), - std::invalid_argument); - } - - SECTION("infinite_exponent equals IEEE 754 positive infinity") { - REQUIRE(infinite_exponent == std::numeric_limits::infinity()); - REQUIRE(std::isinf(infinite_exponent)); - REQUIRE(infinite_exponent > 0.0); - } +TEST_CASE("infinite_exponent sentinel", "[engine][1-body][validation]") { + REQUIRE(libint2::infinite_exponent == + std::numeric_limits::infinity()); + REQUIRE(std::isinf(libint2::infinite_exponent)); + REQUIRE(libint2::infinite_exponent > 0.0); } TEST_CASE("make_q_gau_data factory validation", diff --git a/include/libint2/basis.h.in b/include/libint2/basis.h.in index 8bc2c26ca..075b8af99 100644 --- a/include/libint2/basis.h.in +++ b/include/libint2/basis.h.in @@ -854,7 +854,7 @@ namespace libint2 { // 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) + 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)) diff --git a/include/libint2/shell.h b/include/libint2/shell.h index 96adbeaa3..02ee35704 100644 --- a/include/libint2/shell.h +++ b/include/libint2/shell.h @@ -1354,25 +1354,6 @@ struct GaussianPotentialPrimitive { /// corrections, erf/erfc attenuations, or any combination thereof. using GaussianPotentialData = std::vector; -/// 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: @@ -1385,8 +1366,8 @@ inline void validate_gaussian_potential_primitive( /// 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. -/// When constructing primitives by hand, call -/// validate_gaussian_potential_primitive() on each to catch invalid inputs: +/// 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()); From fc1cdffd47e91a9da3f6cac653d8884d868c9aed Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Fri, 17 Apr 2026 12:41:14 -0400 Subject: [PATCH 9/9] reject zero exponents everywhere --- include/libint2/basis.h.in | 2 +- include/libint2/boys.h | 4 ++-- include/libint2/shell.h | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/libint2/basis.h.in b/include/libint2/basis.h.in index 075b8af99..8bc2c26ca 100644 --- a/include/libint2/basis.h.in +++ b/include/libint2/basis.h.in @@ -854,7 +854,7 @@ namespace libint2 { // 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) + 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)) diff --git a/include/libint2/boys.h b/include/libint2/boys.h index 9a26e8786..fcd11ad7d 100644 --- a/include/libint2/boys.h +++ b/include/libint2/boys.h @@ -2201,8 +2201,8 @@ struct q_gau_gm_eval : private detail::CoreEvalScratch> { 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(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"); diff --git a/include/libint2/shell.h b/include/libint2/shell.h index 02ee35704..d0c52939c 100644 --- a/include/libint2/shell.h +++ b/include/libint2/shell.h @@ -1366,7 +1366,7 @@ using GaussianPotentialData = std::vector; /// 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. -/// Each primitive must satisfy: exponent >= 0 or exponent == +/// Each primitive must satisfy: exponent > 0 or exponent == /// libint2::infinite_exponent; coefficient must be finite (not NaN or +/-inf). /// @code /// using libint2::infinite_exponent;