From 0424882b98e3217a363f0bb7bfa37063a9893ae8 Mon Sep 17 00:00:00 2001 From: Kshitij Surjuse Date: Wed, 25 Mar 2026 11:44:09 -0400 Subject: [PATCH 1/5] 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/5] 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/5] 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/5] 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/5] 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