Skip to content

Adds more error models for primitive screening. #169

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
ryanmrichard wants to merge 3 commits into from
Closed

Adds more error models for primitive screening. #169

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
391a189
adds more error models
ryanmrichard Jul 6, 2026
6512ca7
code cleanup
ryanmrichard Jul 9, 2026
9dc3d4e
undo deleting T
ryanmrichard Jul 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
191 changes: 106 additions & 85 deletions src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <cmath>
#include <integrals/property_types.hpp>
#include <wtf/wtf.hpp>

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<pt>();
description(desc);
// TODO: Add citation for Chemist paper
add_submodule<decontract_pt>("Decontract Basis Set");
add_submodule<eri4_pt>("ERI4");
add_submodule<eri4_pt>("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<decontract_pt>(bra_basis);
const auto& ket_prims = to_prims_mod.run_as<decontract_pt>(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<eri4_pt>(mnls);
chemist::braket::BraKet mnls(bra_pair, v_ee, bra_pair);
const auto& prim4 = submods.at("Raw Primitive ERI4").run_as<eri4_pt>(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<float_type> 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<decltype(n_bra_prims)>; // Type of indices
using index_array = std::array<iter_type, 4>; // Type of a set of 4 indices
using index_vector = std::vector<iter_type>; // Type of a vector of indices
using iter_type = std::size_t;
using index_vector = std::vector<iter_type>;
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<float_type>(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);
}
Expand Down
49 changes: 49 additions & 0 deletions src/integrals/libint/detail_/primitive_pair_estimators.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double,3>` 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<double, 3>;
using vector_t = std::vector<center_t>;
using matrix_t = std::vector<vector_t>;
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_
14 changes: 11 additions & 3 deletions src/integrals/libint/libint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,18 +89,26 @@ 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");
mm.change_submod("Primitive Contractor ERI4", "Raw Primitive ERI4",
"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<LIBINT(bra, op, ket)>(key)
Expand Down
Loading
Loading