From 391a18917c79dae71da562e4d52b19a6b8f81e30 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Mon, 6 Jul 2026 09:57:56 -0500 Subject: [PATCH 1/3] adds more error models --- .../cauchy_schwarz_primitive_estimator.cpp | 191 ++++++++++-------- .../detail_/primitive_pair_estimators.hpp | 49 +++++ src/integrals/libint/libint.cpp | 14 +- .../libint/primitive_error_model.cpp | 137 ++++++++++++- .../libint/primitive_error_model.cpp | 57 ++++++ 5 files changed, 351 insertions(+), 97 deletions(-) diff --git a/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp index 2a6718ea..a1763f2b 100644 --- a/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp +++ b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -14,134 +14,155 @@ * limitations under the License. */ -#include "../utils/rank2_shell_norm.hpp" #include "detail_/make_libint_basis_set.hpp" #include "libint.hpp" #include #include +#include namespace integrals::libint { namespace { -const auto desc = ""; +const auto desc = R"( +CauchySchwarz Primitive Pair Estimator +======================================= + +For each primitive shell pair (pi, pj) from bra_basis x ket_basis, returns: + + cspe[pi][pj] = |c_pi * c_pj| * Q_CS(pi, pj) + +where + - c_pi = d_pi * N_pi / sqrt(shell_norm_i) is the libint-renormalized + contraction coefficient (embed_normalization=true), identical to the + PrimitiveNormalization output used by PrimitiveContractor. + - Q_CS(pi, pj) = sqrt(max_{a,b} (pi_a pj_b | pi_a pj_b)_raw) + is the shell-level Cauchy-Schwarz factor: the square root of the largest + diagonal self-pair ERI over all angular-momentum components a of pi and b of + pj, computed with libint normalization disabled so no N-factor enters. + +By the ERI Cauchy-Schwarz inequality applied to specific AO components +(m_i, m_j): + + |(pi_{m_i} pj_{m_j} | pk_{m_k} pl_{m_l})_raw| + ≤ sqrt((pi_{m_i} pj_{m_j} | pi_{m_i} pj_{m_j})_raw) + * sqrt((pk_{m_k} pl_{m_l} | pk_{m_k} pl_{m_l})_raw) + ≤ Q_CS(pi, pj) * Q_CS(pk, pl) + +so the product cspe[pi][pj] * cspe[pk][pl] is a rigorous upper bound on the +magnitude of any single primitive-AO ERI contribution to a contracted AO element. +)"; } // namespace -using decontract_pt = integrals::property_types::DecontractBasisSet; -using eri4_pt = simde::ERI4; -using pt = integrals::property_types::PrimitivePairEstimator; +using eri4_pt = simde::ERI4; +using pt = integrals::property_types::PrimitivePairEstimator; MODULE_CTOR(CauchySchwarzPrimitiveEstimator) { satisfies_property_type(); description(desc); - // TODO: Add citation for Chemist paper - add_submodule("Decontract Basis Set"); - add_submodule("ERI4"); + add_submodule("Raw Primitive ERI4"); } MODULE_RUN(CauchySchwarzPrimitiveEstimator) { const auto&& [bra_basis, ket_basis] = pt::unwrap_inputs(inputs); - auto& to_prims_mod = submods.at("Decontract Basis Set"); - const auto& bra_prims = to_prims_mod.run_as(bra_basis); - const auto& ket_prims = to_prims_mod.run_as(ket_basis); - const auto n_bra_prims = bra_prims.n_primitives(); - const auto n_ket_prims = ket_prims.n_primitives(); - - // Should always be true, but we check for sanity - assert(n_bra_prims == bra_basis.n_primitives()); - assert(n_ket_prims == ket_basis.n_primitives()); + const auto n_bra_prims = bra_basis.n_primitives(); + const auto n_ket_prims = ket_basis.n_primitives(); - // TODO: We only need the hyper diagonal, so this is very wasteful - simde::type::aos_squared bra(bra_prims, ket_prims); + // Compute the raw (unnormalized) self-pair ERIs for the primitive shells. + // Raw Primitive ERI4 decontracts internally and disables libint + // normalization, so prim4[a,b,a',b'] = raw_ERI without any N-factor. + simde::type::aos bra_aos(bra_basis); + simde::type::aos ket_aos(ket_basis); + simde::type::aos_squared bra_pair(bra_aos, ket_aos); simde::type::v_ee_type v_ee{}; - chemist::braket::BraKet mnls(bra, v_ee, bra); - const auto& prim4 = submods.at("ERI4").run_as(mnls); + chemist::braket::BraKet mnls(bra_pair, v_ee, bra_pair); + const auto& prim4 = submods.at("Raw Primitive ERI4").run_as(mnls); - // TODO: Make our basis set normalize itself. + // Contracted coefficients with libint normalization embedded (default): + // coeff[p] = d_p * N_p / sqrt(contracted_shell_norm) + // This matches the PrimitiveNormalization module used by + // PrimitiveContractor. auto bra_libint = detail_::make_libint_basis_set(bra_basis); auto ket_libint = detail_::make_libint_basis_set(ket_basis); using tensorwrapper::buffer::make_contiguous; const auto& eris = make_contiguous(prim4.buffer()); - // TODO: Use floating point type of the basis sets using float_type = double; std::vector data(n_bra_prims * n_ket_prims, 0.0); tensorwrapper::shape::Smooth shape({n_bra_prims, n_ket_prims}); tensorwrapper::buffer::Contiguous buffer(std::move(data), shape); - using iter_type = std::decay_t; // Type of indices - using index_array = std::array; // Type of a set of 4 indices - using index_vector = std::vector; // Type of a vector of indices + using iter_type = std::size_t; + using index_vector = std::vector; + using wtf::fp::float_cast; - index_array ao_offsets{0, 0, 0, 0}; - index_array naos{0, 0, 0, 0}; - index_vector shell{0, 0}; - index_vector prim{0, 0}; - index_vector prim_offsets{0, 0}; index_vector abs_prim{0, 0}; - for(shell[0] = 0; shell[0] < bra_basis.n_shells(); ++shell[0]) { - const auto& bra_shell = bra_libint.at(shell[0]); - assert(bra_shell.contr.size() == 1); // No general contraction support - const auto& bra_coeff = bra_shell.contr[0].coeff; - const auto n_prims_bra_shell = bra_coeff.size(); - - ao_offsets[0] = 0; - ao_offsets[2] = 0; - for(prim[0] = 0; prim[0] < n_prims_bra_shell; ++prim[0]) { - const auto c_i = std::fabs(bra_coeff[prim[0]]); - abs_prim[0] = prim_offsets[0] + prim[0]; - naos[0] = bra_basis.shell(shell[0]).size(); - naos[2] = naos[0]; - - prim_offsets[1] = 0; - ao_offsets[1] = 0; - ao_offsets[3] = 0; - - for(shell[1] = 0; shell[1] < ket_basis.n_shells(); ++shell[1]) { - const auto& ket_shell = ket_libint.at(shell[1]); - assert(ket_shell.contr.size() == 1); // No general contractions - const auto& ket_coeff = ket_shell.contr[0].coeff; - const auto n_prims_ket_shell = ket_coeff.size(); - - for(prim[1] = 0; prim[1] < n_prims_ket_shell; ++prim[1]) { - const auto c_j = std::fabs(ket_coeff[prim[1]]); - abs_prim[1] = prim_offsets[1] + prim[1]; + // Accumulated absolute AO offsets into the prim4 tensor dimensions. + // For each contracted shell s with n_prims primitives and n_aos AO + // components, primitive p occupies AOs [shell_ao_offset + p*n_aos, + // shell_ao_offset + (p+1)*n_aos). + std::size_t bra_prim_offset = 0; + std::size_t bra_ao_offset = 0; - naos[1] = ket_basis.shell(shell[1]).size(); - naos[3] = naos[1]; - - auto C_ij = c_i * c_j; - - // ao_offsets/Naos needs to respectively be the offset for - // the first "AO" and the number of "AOs" in the - // decontracted ijij shell quartet - - auto shell_norm = - utils::rank2_shell_norm(eris, ao_offsets, naos); - buffer.set_elem(abs_prim, C_ij * shell_norm); - - ao_offsets[1] += naos[1]; - ao_offsets[3] += naos[3]; - - } // loop over ket primitives - - prim_offsets[1] += n_prims_ket_shell; - } // loop over ket shells - - ao_offsets[0] += naos[0]; - ao_offsets[2] += naos[2]; + for(std::size_t s0 = 0; s0 < bra_basis.n_shells(); ++s0) { + const auto& bra_shell = bra_libint.at(s0); + assert(bra_shell.contr.size() == 1); // No general contraction support + const auto& bra_coeff = bra_shell.contr[0].coeff; + const auto n_prims_s0 = bra_coeff.size(); + const auto naos_s0 = bra_basis.shell(s0).size(); - } // loop over bra primitives + for(std::size_t p0 = 0; p0 < n_prims_s0; ++p0) { + const auto c_i = std::fabs(bra_coeff[p0]); + abs_prim[0] = bra_prim_offset + p0; + const auto off_pi = bra_ao_offset + p0 * naos_s0; - prim_offsets[0] += n_prims_bra_shell; + std::size_t ket_prim_offset = 0; + std::size_t ket_ao_offset = 0; - } // loop over bra shells + for(std::size_t s1 = 0; s1 < ket_basis.n_shells(); ++s1) { + const auto& ket_shell = ket_libint.at(s1); + assert(ket_shell.contr.size() == 1); // No general contractions + const auto& ket_coeff = ket_shell.contr[0].coeff; + const auto n_prims_s1 = ket_coeff.size(); + const auto naos_s1 = ket_basis.shell(s1).size(); + + for(std::size_t p1 = 0; p1 < n_prims_s1; ++p1) { + const auto c_j = std::fabs(ket_coeff[p1]); + abs_prim[1] = ket_prim_offset + p1; + const auto off_pj = ket_ao_offset + p1 * naos_s1; + + // Q_CS(pi,pj) = sqrt(max_{a,b} (pi_a pj_b | pi_a pj_b)_raw) + // This is the correct Cauchy-Schwarz factor: by the ERI + // positivity, (pi_a pj_b | pi_a pj_b) >= 0, and by C-S: + // |(pi_a pj_b | pk_c pl_d)| <= sqrt((pi_a pj_b | pi_a + // pj_b)) + // * sqrt((pk_c pl_d | pk_c + // pl_d)) + // <= Q_CS(pi,pj) * Q_CS(pk,pl). + double max_diag = 0.0; + for(iter_type a = 0; a < naos_s0; ++a) { + for(iter_type b = 0; b < naos_s1; ++b) { + index_vector idx4 = {off_pi + a, off_pj + b, + off_pi + a, off_pj + b}; + const auto val = + float_cast(eris.get_elem(idx4)); + max_diag = std::max(max_diag, val); + } + } + buffer.set_elem(abs_prim, c_i * c_j * std::sqrt(max_diag)); + } + ket_prim_offset += n_prims_s1; + ket_ao_offset += n_prims_s1 * naos_s1; + } + } + bra_prim_offset += n_prims_s0; + bra_ao_offset += n_prims_s0 * naos_s0; + } simde::type::tensor rv(shape, std::move(buffer)); - auto result = results(); return pt::wrap_results(result, rv); } diff --git a/src/integrals/libint/detail_/primitive_pair_estimators.hpp b/src/integrals/libint/detail_/primitive_pair_estimators.hpp index 08604657..2ee845f3 100644 --- a/src/integrals/libint/detail_/primitive_pair_estimators.hpp +++ b/src/integrals/libint/detail_/primitive_pair_estimators.hpp @@ -166,4 +166,53 @@ inline auto fine_k_ij(const simde::type::ao_basis_set& basis0, } return K; } +/** @brief Gaussian product centers + * @f$\mathbf{P}_{ij} = (\alpha_i \mathbf{A} + \alpha_j \mathbf{B}) / + * \gamma_{ij}@f$ for each primitive pair. + * + * @param[in] basis0 First basis (rows). + * @param[in] basis1 Second basis (columns). + * @return Matrix of shape `n_prims0` x `n_prims1`, each entry an + * `std::array` holding (x, y, z) of the product center. + */ +inline auto product_centers_ij(const simde::type::ao_basis_set& basis0, + const simde::type::ao_basis_set& basis1) { + auto nprims0 = basis0.n_primitives(); + auto nprims1 = basis1.n_primitives(); + auto gamma = gamma_ij(basis0, basis1); + using center_t = std::array; + using vector_t = std::vector; + using matrix_t = std::vector; + matrix_t P(nprims0, vector_t(nprims1, {0.0, 0.0, 0.0})); + for(std::size_t i = 0; i < nprims0; ++i) { + auto alpha0 = basis0.primitive(i).exponent(); + auto r0 = basis0.primitive(i).center(); + for(std::size_t j = 0; j < nprims1; ++j) { + auto alpha1 = basis1.primitive(j).exponent(); + auto r1 = basis1.primitive(j).center(); + auto gij = gamma[i][j]; + P[i][j] = {(alpha0 * r0.x() + alpha1 * r1.x()) / gij, + (alpha0 * r0.y() + alpha1 * r1.y()) / gij, + (alpha0 * r0.z() + alpha1 * r1.z()) / gij}; + } + } + return P; +} + +/** @brief Upper bound on the @f$F_0(T)@f$ Boys function. + * + * Uses @f$F_0(T) \leq \min\!\left(1,\, + * \tfrac{\sqrt{\pi}}{2\sqrt{T}}\right)@f$, which is exact in the large-@f$T@f$ + * limit and the trivial bound 1 for + * @f$T \leq \pi/4@f$. + * + * @param[in] T Non-negative argument @f$T = \frac{\gamma_{ij}\gamma_{kl}} + * {\gamma_{ij}+\gamma_{kl}} R_{PQ}^2@f$. + * @return Upper bound on @f$F_0(T)@f$ in [0, 1]. + */ +inline double boys_f0_upper_bound(double T) { + if(T <= 0.0) return 1.0; + return std::min(1.0, 0.5 * std::sqrt(M_PI / T)); +} + } // namespace integrals::libint::detail_ diff --git a/src/integrals/libint/libint.cpp b/src/integrals/libint/libint.cpp index f453ddfb..1acce970 100644 --- a/src/integrals/libint/libint.cpp +++ b/src/integrals/libint/libint.cpp @@ -89,11 +89,17 @@ EXTERN_LIBINT(aos_squared, v_ee_type, aos_squared); #undef EXTERN_LIBINT void set_defaults(pluginplay::ModuleManager& mm) { - mm.change_submod("CauchySchwarz Estimator", "Decontract Basis Set", - "Decontract Basis Set"); mm.copy_module("ERI4", "Benchmark ERI4"); mm.change_input("Benchmark ERI4", "Threshold", 1.0E-16); - mm.change_submod("CauchySchwarz Estimator", "ERI4", "Benchmark ERI4"); + // Separate Raw Primitive ERI4 copy for Schwarz self-pair computation. + // Threshold=0 ensures no self-pair ERI is rounded to zero, which would + // cause a false Q_CS=0 and violate the bound for tiny-but-nonzero ERIs. + mm.copy_module("Raw Primitive ERI4", "Schwarz Raw Primitive ERI4"); + mm.change_input("Schwarz Raw Primitive ERI4", "Threshold", 0.0); + mm.change_submod("Schwarz Raw Primitive ERI4", "Decontract Basis Set", + "Decontract Basis Set"); + mm.change_submod("CauchySchwarz Estimator", "Raw Primitive ERI4", + "Schwarz Raw Primitive ERI4"); mm.change_submod("Analytic Error", "ERI4s", "Benchmark ERI4"); mm.change_submod("Raw Primitive ERI4", "Decontract Basis Set", "Decontract Basis Set"); @@ -101,6 +107,8 @@ void set_defaults(pluginplay::ModuleManager& mm) { "Raw Primitive ERI4"); mm.change_submod("Primitive Contractor ERI4", "Primitive Normalization", "Primitive Normalization"); + mm.change_submod("Primitive Error Model", "CauchySchwarz Estimator", + "CauchySchwarz Estimator"); } #define LOAD_LIBINT(bra, op, ket, key) mm.add_module(key) diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp index 48255480..76f8cf29 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -19,8 +19,10 @@ #include "libint.hpp" #include #include +#include #include #include +#include namespace integrals::libint { namespace { @@ -35,9 +37,21 @@ same coarse / fine gates as libint `ScreeningMethod::Original` (coarse `coarse_k_ij`, fine `fine_k_ij` with `gamma_ij`). For each decontracted quartet that would be skipped by the contractor, this -module accumulates one of three per-quartet estimates into the corresponding -contracted AO element: fixed `Tolerance`, coarse pair product -`K_ij * K_kl`, or fine metric `|Q_ij Q_kl| / sqrt(gamma_ij + gamma_kl)`. +module accumulates one of seven per-quartet estimates into the corresponding +contracted AO element: fixed `Tolerance`, coarse pair product `K_ij * K_kl`, +fine metric `|Q_ij Q_kl| / sqrt(gamma_ij + gamma_kl)`, `FineBoys` which +multiplies the fine metric by an upper bound on F_0(T) (Boys-function +diagnostic, NOT a rigorous upper bound for higher angular momenta), `Schwarz` +which uses sqrt(||(ij|ij)||_F) * sqrt(||(kl|kl)||_F) from the CauchySchwarz +submodule — a rigorous upper bound that correctly captures all angular-momentum +effects, `SchwarzBoys` which multiplies the Schwarz product by the same F_0(T) +upper bound — tighter than plain Schwarz for well-separated charge distributions +(rigorous for s-only quartets; empirically validated for higher angular +momenta), or `SchwarzGF` which additionally multiplies by the exponent-mismatch +factor G(p,q) = (4 p q / (p+q)^2)^(1/4) with p = gamma_ij, q = gamma_kl. For an +s-only quartet Schwarz * G(p,q) * F_0(T) is the EXACT integral (zero +overestimation); for higher angular momenta it is a non-rigorous tightened +estimate. )"; /** @return True iff `PrimitiveContractor` would `continue` (skip) this quartet. @@ -54,25 +68,61 @@ inline bool primitive_quartet_skipped(double K_ij, double K_kl, double Q_ij, return false; } -enum class ErrorEstimateKind { Tolerance, Coarse, Fine }; +enum class ErrorEstimateKind { + Tolerance, + Coarse, + Fine, + FineBoys, + Schwarz, + SchwarzBoys, + SchwarzGF +}; inline ErrorEstimateKind parse_error_estimate(const std::string& s) { if(s == "Tolerance") return ErrorEstimateKind::Tolerance; if(s == "Coarse") return ErrorEstimateKind::Coarse; if(s == "Fine") return ErrorEstimateKind::Fine; + if(s == "FineBoys") return ErrorEstimateKind::FineBoys; + if(s == "Schwarz") return ErrorEstimateKind::Schwarz; + if(s == "SchwarzBoys") return ErrorEstimateKind::SchwarzBoys; + if(s == "SchwarzGF") return ErrorEstimateKind::SchwarzGF; throw std::invalid_argument( "Primitive Error Model: \"Error estimate\" must be \"Tolerance\", " - "\"Coarse\", or \"Fine\""); + "\"Coarse\", \"Fine\", \"FineBoys\", \"Schwarz\", \"SchwarzBoys\", or " + "\"SchwarzGF\""); +} + +/** @brief Exponent-mismatch factor + * @f$G(p,q) = \left(4pq / (p+q)^2\right)^{1/4} = + * \sqrt{2\sqrt{pq}/(p+q)} + * @f$. + * + * For an s-type primitive quartet the exact ratio of the true ERI to the + * Cauchy-Schwarz product is @f$G(p,q)\,F_0(T)@f$, where @f$p=\gamma_{ij}@f$ + * and @f$q=\gamma_{kl}@f$. By AM-GM @f$G \le 1@f$, with equality iff + * @f$p=q@f$; it measures the tight-vs-diffuse mismatch between the bra and ket + * pairs that the plain Schwarz bound ignores. + */ +inline double gf_exponent_factor(double gamma_ij, double gamma_kl) { + const double s = gamma_ij + gamma_kl; + return std::pow(4.0 * gamma_ij * gamma_kl / (s * s), 0.25); } inline double skip_increment(ErrorEstimateKind kind, double thresh, double K_ij, double K_kl, double Q_ij, double Q_kl, - double gamma_ij, double gamma_kl) { + double gamma_ij, double gamma_kl, double T = 0.0) { switch(kind) { case ErrorEstimateKind::Tolerance: return thresh; case ErrorEstimateKind::Coarse: return K_ij * K_kl; case ErrorEstimateKind::Fine: return std::abs(Q_ij * Q_kl / std::sqrt(gamma_ij + gamma_kl)); + case ErrorEstimateKind::FineBoys: + return std::abs(Q_ij * Q_kl / std::sqrt(gamma_ij + gamma_kl)) * + detail_::boys_f0_upper_bound(T); + case ErrorEstimateKind::Schwarz: + case ErrorEstimateKind::SchwarzBoys: + case ErrorEstimateKind::SchwarzGF: + return 0.0; // handled separately in the quartet loop } return 0.0; } @@ -80,6 +130,7 @@ inline double skip_increment(ErrorEstimateKind kind, double thresh, double K_ij, } // namespace using eri4_pt = simde::ERI4; +using ppt = integrals::property_types::PrimitivePairEstimator; using pt = integrals::property_types::Uncertainty; MODULE_CTOR(PrimitiveErrorModel) { @@ -91,7 +142,20 @@ MODULE_CTOR(PrimitiveErrorModel) { .set_description( "Per skipped primitive quartet: \"Tolerance\" adds the screening " "threshold; \"Coarse\" adds K_ij*K_kl; \"Fine\" adds the fine-screen " - "metric |Q_ij Q_kl|/sqrt(gamma_ij+gamma_kl)."); + "metric |Q_ij Q_kl|/sqrt(gamma_ij+gamma_kl); \"FineBoys\" multiplies " + "the Fine metric by an upper bound on F_0(T) (diagnostic only); " + "\"Schwarz\" uses the Cauchy-Schwarz bound sqrt(||(ij|ij)||_F) * " + "sqrt(||(kl|kl)||_F) — rigorous for all angular momenta; " + "\"SchwarzBoys\" multiplies the Schwarz product by an upper bound on " + "F_0(T) — tighter for well-separated charge distributions (rigorous " + "for s-only quartets, empirically validated for higher angular " + "momenta); \"SchwarzGF\" additionally multiplies by the " + "exponent-mismatch factor G(p,q)=(4pq/(p+q)^2)^(1/4) — exact for " + "s-only " + "quartets, non-rigorous tightened estimate for higher angular " + "momenta."); + + add_submodule("CauchySchwarz Estimator"); } MODULE_RUN(PrimitiveErrorModel) { @@ -114,6 +178,35 @@ MODULE_RUN(PrimitiveErrorModel) { const auto Q_bra = detail_::fine_k_ij(bs0, bs1); const auto Q_ket = detail_::fine_k_ij(bs2, bs3); + const bool need_boys = (kind == ErrorEstimateKind::FineBoys || + kind == ErrorEstimateKind::SchwarzBoys || + kind == ErrorEstimateKind::SchwarzGF); + const bool need_schwarz = (kind == ErrorEstimateKind::Schwarz || + kind == ErrorEstimateKind::SchwarzBoys || + kind == ErrorEstimateKind::SchwarzGF); + + const auto P_bra = need_boys ? + detail_::product_centers_ij(bs0, bs1) : + decltype(detail_::product_centers_ij(bs0, bs1)){}; + const auto P_ket = need_boys ? + detail_::product_centers_ij(bs2, bs3) : + decltype(detail_::product_centers_ij(bs2, bs3)){}; + + // Schwarz tensors must outlive the spans derived from them. + simde::type::tensor S_bra_tensor, S_ket_tensor; + std::span schwarz_bra_data, schwarz_ket_data; + const std::size_t n1_prims = bs1.n_primitives(); + const std::size_t n3_prims = bs3.n_primitives(); + if(need_schwarz) { + auto& cs_mod = submods.at("CauchySchwarz Estimator"); + S_bra_tensor = cs_mod.run_as(bs0, bs1); + S_ket_tensor = cs_mod.run_as(bs2, bs3); + schwarz_bra_data = + tensorwrapper::buffer::get_raw_data(S_bra_tensor.buffer()); + schwarz_ket_data = + tensorwrapper::buffer::get_raw_data(S_ket_tensor.buffer()); + } + auto map0 = utils::build_prim_ao_to_cgto_map(bs0); auto map1 = utils::build_prim_ao_to_cgto_map(bs1); auto map2 = utils::build_prim_ao_to_cgto_map(bs2); @@ -169,8 +262,34 @@ MODULE_RUN(PrimitiveErrorModel) { continue; } - const double inc = skip_increment( - kind, tol, K_ij, K_kl, Q_ij, Q_kl, gamma_ij, gamma_kl); + double T = 0.0; + if(need_boys) { + const auto& Pij = P_bra[pi][pj]; + const auto& Pkl = P_ket[pk][pl]; + const double dx = Pij[0] - Pkl[0]; + const double dy = Pij[1] - Pkl[1]; + const double dz = Pij[2] - Pkl[2]; + T = gamma_ij * gamma_kl / (gamma_ij + gamma_kl) * + (dx * dx + dy * dy + dz * dz); + } + + double inc = 0.0; + if(kind == ErrorEstimateKind::SchwarzBoys) { + inc = schwarz_bra_data[pi * n1_prims + pj] * + schwarz_ket_data[pk * n3_prims + pl] * + detail_::boys_f0_upper_bound(T); + } else if(kind == ErrorEstimateKind::SchwarzGF) { + inc = schwarz_bra_data[pi * n1_prims + pj] * + schwarz_ket_data[pk * n3_prims + pl] * + gf_exponent_factor(gamma_ij, gamma_kl) * + detail_::boys_f0_upper_bound(T); + } else if(need_schwarz) { + inc = schwarz_bra_data[pi * n1_prims + pj] * + schwarz_ket_data[pk * n3_prims + pl]; + } else { + inc = skip_increment(kind, tol, K_ij, K_kl, Q_ij, Q_kl, + gamma_ij, gamma_kl, T); + } raw_buffer[ao_offset(mu, nu, lam, sig)] += inc; } } diff --git a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp index 08cf56cc..6a86483a 100644 --- a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -29,9 +29,11 @@ using namespace integrals::testing; namespace { +using integrals::libint::detail_::boys_f0_upper_bound; using integrals::libint::detail_::coarse_k_ij; using integrals::libint::detail_::fine_k_ij; using integrals::libint::detail_::gamma_ij; +using integrals::libint::detail_::product_centers_ij; bool quartet_skipped(double K_ij, double K_kl, double Q_ij, double Q_kl, double gamma_ij, double gamma_kl, double thresh) { @@ -146,6 +148,61 @@ TEST_CASE("PrimitiveErrorModel") { Catch::Approx(ref.sum_fine).margin(1e-10)); } + SECTION("SchwarzBoys ≤ Schwarz element-wise (H2 STO-3G)") { + auto aobs = h2_sto3g_basis_set(); + simde::type::aos_squared bra(aobs, aobs); + simde::type::aos_squared ket(aobs, aobs); + chemist::braket::BraKet mnls(bra, v_ee, ket); + + const double tol = 1e-12; + auto t_schwarz = run_mode(mnls, tol, "Schwarz"); + auto t_schwarzboys = run_mode(mnls, tol, "SchwarzBoys"); + + using tensorwrapper::buffer::get_raw_data; + const auto buf_s = get_raw_data(t_schwarz.buffer()); + const auto buf_sb = get_raw_data(t_schwarzboys.buffer()); + + REQUIRE(buf_s.size() == buf_sb.size()); + bool any_strictly_smaller = false; + for(std::size_t e = 0; e < buf_s.size(); ++e) { + REQUIRE(buf_sb[e] <= buf_s[e] + 1e-14); + if(buf_sb[e] < buf_s[e] - 1e-14) any_strictly_smaller = true; + } + // For H2 there are separated primitive pairs, so F0 < 1 for some. + REQUIRE(any_strictly_smaller); + REQUIRE(buffer_sum(t_schwarzboys) < buffer_sum(t_schwarz)); + } + + SECTION("SchwarzGF ≤ SchwarzBoys ≤ Schwarz element-wise (H2 STO-3G)") { + auto aobs = h2_sto3g_basis_set(); + simde::type::aos_squared bra(aobs, aobs); + simde::type::aos_squared ket(aobs, aobs); + chemist::braket::BraKet mnls(bra, v_ee, ket); + + const double tol = 1e-12; + auto t_schwarz = run_mode(mnls, tol, "Schwarz"); + auto t_schwarzboys = run_mode(mnls, tol, "SchwarzBoys"); + auto t_schwarzgf = run_mode(mnls, tol, "SchwarzGF"); + + using tensorwrapper::buffer::get_raw_data; + const auto buf_s = get_raw_data(t_schwarz.buffer()); + const auto buf_sb = get_raw_data(t_schwarzboys.buffer()); + const auto buf_gf = get_raw_data(t_schwarzgf.buffer()); + + REQUIRE(buf_s.size() == buf_gf.size()); + bool any_strictly_smaller = false; + for(std::size_t e = 0; e < buf_s.size(); ++e) { + // G(p,q) ≤ 1 so SchwarzGF ≤ SchwarzBoys ≤ Schwarz. + REQUIRE(buf_gf[e] <= buf_sb[e] + 1e-14); + REQUIRE(buf_gf[e] <= buf_s[e] + 1e-14); + if(buf_gf[e] < buf_sb[e] - 1e-14) any_strictly_smaller = true; + } + // H2/STO-3G mixes exponents within a contraction, so some pairs have + // p != q and G(p,q) < 1. + REQUIRE(any_strictly_smaller); + REQUIRE(buffer_sum(t_schwarzgf) < buffer_sum(t_schwarzboys)); + } + SECTION("invalid Error estimate throws") { auto aobs = h2_sto3g_basis_set(); simde::type::aos_squared bra(aobs, aobs); From 6512ca76d10802f8f8c8ecf0a48ce1119955b9e3 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Thu, 9 Jul 2026 11:38:17 -0500 Subject: [PATCH 2/3] code cleanup --- .../libint/primitive_error_model.cpp | 27 ++----------------- 1 file changed, 2 insertions(+), 25 deletions(-) diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp index 76f8cf29..97e1c566 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -42,7 +42,7 @@ contracted AO element: fixed `Tolerance`, coarse pair product `K_ij * K_kl`, fine metric `|Q_ij Q_kl| / sqrt(gamma_ij + gamma_kl)`, `FineBoys` which multiplies the fine metric by an upper bound on F_0(T) (Boys-function diagnostic, NOT a rigorous upper bound for higher angular momenta), `Schwarz` -which uses sqrt(||(ij|ij)||_F) * sqrt(||(kl|kl)||_F) from the CauchySchwarz +which uses sqrt(||(ij|ij)||) * sqrt(||(kl|kl)||) from the CauchySchwarz submodule — a rigorous upper bound that correctly captures all angular-momentum effects, `SchwarzBoys` which multiplies the Schwarz product by the same F_0(T) upper bound — tighter than plain Schwarz for well-separated charge distributions @@ -139,21 +139,7 @@ MODULE_CTOR(PrimitiveErrorModel) { add_input("Error estimate") .set_default("Tolerance") - .set_description( - "Per skipped primitive quartet: \"Tolerance\" adds the screening " - "threshold; \"Coarse\" adds K_ij*K_kl; \"Fine\" adds the fine-screen " - "metric |Q_ij Q_kl|/sqrt(gamma_ij+gamma_kl); \"FineBoys\" multiplies " - "the Fine metric by an upper bound on F_0(T) (diagnostic only); " - "\"Schwarz\" uses the Cauchy-Schwarz bound sqrt(||(ij|ij)||_F) * " - "sqrt(||(kl|kl)||_F) — rigorous for all angular momenta; " - "\"SchwarzBoys\" multiplies the Schwarz product by an upper bound on " - "F_0(T) — tighter for well-separated charge distributions (rigorous " - "for s-only quartets, empirically validated for higher angular " - "momenta); \"SchwarzGF\" additionally multiplies by the " - "exponent-mismatch factor G(p,q)=(4pq/(p+q)^2)^(1/4) — exact for " - "s-only " - "quartets, non-rigorous tightened estimate for higher angular " - "momenta."); + .set_description(desc); add_submodule("CauchySchwarz Estimator"); } @@ -263,15 +249,6 @@ MODULE_RUN(PrimitiveErrorModel) { } double T = 0.0; - if(need_boys) { - const auto& Pij = P_bra[pi][pj]; - const auto& Pkl = P_ket[pk][pl]; - const double dx = Pij[0] - Pkl[0]; - const double dy = Pij[1] - Pkl[1]; - const double dz = Pij[2] - Pkl[2]; - T = gamma_ij * gamma_kl / (gamma_ij + gamma_kl) * - (dx * dx + dy * dy + dz * dz); - } double inc = 0.0; if(kind == ErrorEstimateKind::SchwarzBoys) { From 9dc3d4e32f8a6a24a0ee0c5bfb20806521394c0a Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Wed, 15 Jul 2026 08:47:28 -0500 Subject: [PATCH 3/3] undo deleting T --- src/integrals/libint/primitive_error_model.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp index 97e1c566..1ad5cd33 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -249,6 +249,15 @@ MODULE_RUN(PrimitiveErrorModel) { } double T = 0.0; + if(need_boys) { + const auto& Pij = P_bra[pi][pj]; + const auto& Pkl = P_ket[pk][pl]; + const double dx = Pij[0] - Pkl[0]; + const double dy = Pij[1] - Pkl[1]; + const double dz = Pij[2] - Pkl[2]; + T = gamma_ij * gamma_kl / (gamma_ij + gamma_kl) * + (dx * dx + dy * dy + dz * dz); + } double inc = 0.0; if(kind == ErrorEstimateKind::SchwarzBoys) {