diff --git a/src/integrals/libint/detail_/primitive_pair_estimators.hpp b/src/integrals/libint/detail_/primitive_pair_estimators.hpp new file mode 100644 index 00000000..08604657 --- /dev/null +++ b/src/integrals/libint/detail_/primitive_pair_estimators.hpp @@ -0,0 +1,169 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +/** @file primitive_pair_estimators.hpp + * @brief Dense primitive-pair matrices for ERI screening and prefactors. + * + * These helpers build `n_primitives(basis0) x n_primitives(basis1)` matrices + * for two `simde::type::ao_basis_set` arguments. + * + * **Indexing:** Row `i` is the `i`-th primitive in @p basis0 (Chemist + * flattening: shells in declaration order, primitives fastest within each + * shell). Column `j` is the `j`-th primitive in @p basis1. This matches the + * primitive order used by `BlackBoxPrimitiveEstimator` and libint + * `make_libint_basis_set`. + * + * The **coarse** path uses libint-normalized coefficients from + * `make_libint_basis_set` and assumes **one contraction per shell** + * (`contr[0]`), matching the current implementation. + */ + +#pragma once +#include "make_libint_basis_set.hpp" +#include +#include +#include + +namespace integrals::libint::detail_ { + +/** @brief Sum of exponents @f$\gamma_{ij} = \alpha_i + \alpha_j@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`. + */ +inline auto gamma_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(); + using vector_t = std::vector; + using matrix_t = std::vector; + matrix_t gamma(nprims0, vector_t(nprims1, 0.0)); + for(std::size_t i = 0; i < nprims0; ++i) { + auto alpha0 = basis0.primitive(i).exponent(); + for(std::size_t j = 0; j < nprims1; ++j) { + auto alpha1 = basis1.primitive(j).exponent(); + gamma[i][j] = alpha0 + alpha1; + } + } + return gamma; +} + +/** @brief Overlap-style factor @f$\exp(-\rho_{ij} R_{AB}^2)@f$ with + * @f$\rho_{ij} = \alpha_i \alpha_j / \gamma_{ij}@f$. Primitive centers + * come from Chemist geometry. + * + * @param[in] basis0 First basis (rows). + * @param[in] basis1 Second basis (columns). + * @return Matrix of shape `basis0.n_primitives()` x `basis1.n_primitives()`. + */ +inline auto k_ij(const simde::type::ao_basis_set& basis0, + const simde::type::ao_basis_set& basis1) { + auto distance_squared = [](auto&& a, auto&& b) { + auto dx = a.x() - b.x(); + auto dy = a.y() - b.y(); + auto dz = a.z() - b.z(); + return dx * dx + dy * dy + dz * dz; + }; + + auto nprims0 = basis0.n_primitives(); + auto nprims1 = basis1.n_primitives(); + using vector_t = std::vector; + using matrix_t = std::vector; + auto gamma = gamma_ij(basis0, basis1); + matrix_t K(nprims0, vector_t(nprims1, 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 gamma_ij = gamma[i][j]; + auto rho_ij = alpha0 * alpha1 / gamma_ij; + auto dr2 = distance_squared(r0, r1); + K[i][j] = std::exp(-rho_ij * dr2); + } + } + return K; +} + +/** @brief Coarse screening-style estimate: `k_ij` multiplied by + * @f$|c_i|\,|c_j|@f$ from libint shell coefficients. + * + * Libint uses this pair estimate for "coarse" screening. Compared to the fine + * estimate it neglects the geometric factor. Assumes one contraction per + * shell. + * + * @param[in] basis0 First basis (rows). + * @param[in] basis1 Second basis (columns). + * @return Matrix of shape `basis0.n_primitives()` x `basis1.n_primitives()`. + */ +inline auto coarse_k_ij(const simde::type::ao_basis_set& basis0, + const simde::type::ao_basis_set& basis1) { + auto K = k_ij(basis0, basis1); + auto bs0_libint = make_libint_basis_set(basis0); + auto bs1_libint = make_libint_basis_set(basis1); + std::size_t prim0_offset = 0; + for(const auto& shell0 : bs0_libint) { + auto nprims0 = shell0.nprim(); + for(std::size_t prim0 = 0; prim0 < nprims0; ++prim0) { + auto i = prim0_offset + prim0; + auto coefi = std::fabs(shell0.contr[0].coeff[prim0]); + std::size_t prim1_offset = 0; + for(const auto& shell1 : bs1_libint) { + auto nprims1 = shell1.nprim(); + for(std::size_t prim1 = 0; prim1 < nprims1; ++prim1) { + auto j = prim1_offset + prim1; + auto coefj = std::fabs(shell1.contr[0].coeff[prim1]); + K[i][j] *= coefi * coefj; + } + prim1_offset += nprims1; + } + } + prim0_offset += nprims0; + } + return K; +} + +/** @brief Fine screening-style estimate: `coarse_k_ij` scaled by + * @f$\sqrt{2}\,\pi^{5/4} / \gamma_{ij}@f$. Up to floating-point detail, + * this matches libint `ShellPair` `pp.K` times the same coefficient + * product used in `coarse_k_ij`. + * + * Assumes one contraction per shell. + * + * @note Uses `M_PI` from ``. On MSVC, define `_USE_MATH_DEFINES` before + * including ``, or replace with a portable @f$\pi@f$ if needed. + * + * @param[in] basis0 First basis (rows). + * @param[in] basis1 Second basis (columns). + * @return Matrix of shape `basis0.n_primitives()` x `basis1.n_primitives()`. + */ +inline auto fine_k_ij(const simde::type::ao_basis_set& basis0, + const simde::type::ao_basis_set& basis1) { + auto K = coarse_k_ij(basis0, basis1); + auto gamma = gamma_ij(basis0, basis1); + double prefactor = std::sqrt(2.0) * std::pow(M_PI, 1.25); + for(std::size_t i = 0; i < K.size(); ++i) { + for(std::size_t j = 0; j < K[i].size(); ++j) { + K[i][j] *= prefactor / gamma[i][j]; + } + } + return K; +} +} // namespace integrals::libint::detail_ diff --git a/src/integrals/libint/primitive_contractor.cpp b/src/integrals/libint/primitive_contractor.cpp index 3467c105..d38fd199 100644 --- a/src/integrals/libint/primitive_contractor.cpp +++ b/src/integrals/libint/primitive_contractor.cpp @@ -14,8 +14,14 @@ * limitations under the License. */ +#include "../utils/primitive_index_helpers.hpp" +#include "detail_/make_libint_basis_set.hpp" +#include "detail_/primitive_pair_estimators.hpp" #include "libint.hpp" +#include #include +#include +#include #include namespace integrals::libint { @@ -46,41 +52,20 @@ The algorithm proceeds in three steps: the decontracted basis (primitive index varying fastest within each (shell, ao_component) block). -3. The contracted AO integrals are formed by: +3. Matrices of pair estimates are computed for bra and ket basis sets. K and + K' are respectively coarse estimates for bra and ket basis sets. Q and Q' are + respectively fine estimates for bra and ket basis sets. - AO[mu, nu, lam, sig] = - sum_{i in prims(mu)} sum_{j in prims(nu)} - sum_{k in prims(lam)} sum_{l in prims(sig)} - c0[i] * c1[j] * c2[k] * c3[l] * prim_ERI[i, j, k, l] +4. The contracted AO integrals are formed by summing + c0[i]*c1[j]*c2[k]*c3[l]*prim_ERI[i,j,k,l] with screening: - where prims(mu) is the set of decontracted AO indices that belong to the - same contracted AO mu (same shell, same AO component, all primitives). + Coarse (from `coarse_k_ij`): + - Skip bra pair (i,j) if K[i,j] < t + - Skip ket pair (k,l) if K'[k,l] < t + - Skip quartet if K[i,j] * K'[k,l] <= t (libint keeps only sums > ln(t)) + - Skip if |Q[i,j] * Q'[k,l] / sqrt(gamma_bra + gamma_ket)| < t. )"; -// Build a map from decontracted AO index -> contracted AO index. -// For contracted shell s with n_prims primitives and n_aos AO components, -// the decontracted indices (prim p, ao_component a) all map to the contracted -// AO index first_ao[s] + a. -// NOTE: bs.size() returns the number of centers, not shells. -// bs.n_shells() returns the total number of shells (flattened). -std::vector build_prim_to_ao_map( - const simde::type::ao_basis_set& bs) { - std::vector map; - std::size_t contracted_ao = 0; - for(std::size_t s = 0; s < bs.n_shells(); ++s) { - const auto& shell = bs.shell(s); - const auto n_prims = shell.n_primitives(); - const auto n_aos = shell.size(); - for(std::size_t p = 0; p < n_prims; ++p) { - for(std::size_t a = 0; a < n_aos; ++a) { - map.push_back(contracted_ao + a); - } - } - contracted_ao += n_aos; - } - return map; -} - } // namespace using my_pt = simde::ERI4; @@ -91,10 +76,17 @@ MODULE_CTOR(PrimitiveContractor) { description(desc); add_submodule("Raw Primitive ERI4"); add_submodule("Primitive Normalization"); + add_input("Screening Threshold") + .set_default(1e-16) + .set_description( + "Coarse screening uses the PrimitivePairEstimator matrices; fine " + "screening uses libint ShellPair geometry factors. Matches libint " + "ScreeningMethod::Original."); } MODULE_RUN(PrimitiveContractor) { const auto& [braket] = my_pt::unwrap_inputs(inputs); + const double thresh = inputs.at("Screening Threshold").value(); auto bra = braket.bra(); auto ket = braket.ket(); @@ -115,57 +107,93 @@ MODULE_RUN(PrimitiveContractor) { const auto& c2 = norm_mod.run_as(bs2); const auto& c3 = norm_mod.run_as(bs3); - // Build primitive -> contracted AO index maps - auto map0 = build_prim_to_ao_map(bs0); - auto map1 = build_prim_to_ao_map(bs1); - auto map2 = build_prim_to_ao_map(bs2); - auto map3 = build_prim_to_ao_map(bs3); + // Step 3: coarse screening matrices from PrimitivePairEstimator. + const auto K_bra = detail_::coarse_k_ij(bs0, bs1); + const auto k_ket = detail_::coarse_k_ij(bs2, bs3); + + 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); + auto map3 = utils::build_prim_ao_to_cgto_map(bs3); + + auto pmap0 = utils::build_prim_ao_to_prim_shell_map(bs0); + auto pmap1 = utils::build_prim_ao_to_prim_shell_map(bs1); + auto pmap2 = utils::build_prim_ao_to_prim_shell_map(bs2); + auto pmap3 = utils::build_prim_ao_to_prim_shell_map(bs3); - const Eigen::Index n0 = bs0.n_aos(); - const Eigen::Index n1 = bs1.n_aos(); - const Eigen::Index n2 = bs2.n_aos(); - const Eigen::Index n3 = bs3.n_aos(); + // Number of AOs per basis set + std::array naos{bs0.n_aos(), bs1.n_aos(), bs2.n_aos(), bs3.n_aos()}; - const Eigen::Index np0 = c0.size(); - const Eigen::Index np1 = c1.size(); - const Eigen::Index np2 = c2.size(); - const Eigen::Index np3 = c3.size(); + // Number of primitive shells per basis set + std::array nprims{c0.size(), c1.size(), c2.size(), c3.size()}; - // Get raw data from the primitive tensor via Eigen TensorMap using namespace tensorwrapper; - const auto pdata = buffer::get_raw_data(prim_T.buffer()); - - using prim_map_t = - Eigen::TensorMap>; - prim_map_t prim(pdata.data(), np0, np1, np2, np3); - - // Step 3: contract into AO basis - std::vector ao_data(n0 * n1 * n2 * n3, 0.0); - using ao_map_t = - Eigen::TensorMap>; - ao_map_t ao(ao_data.data(), n0, n1, n2, n3); - - for(Eigen::Index i = 0; i < np0; ++i) { - const double ci = c0[i]; - const Eigen::Index mu = map0[i]; - for(Eigen::Index j = 0; j < np1; ++j) { - const double cij = ci * c1[j]; - const Eigen::Index nu = map1[j]; - for(Eigen::Index k = 0; k < np2; ++k) { - const double cijk = cij * c2[k]; - const Eigen::Index lam = map2[k]; - for(Eigen::Index l = 0; l < np3; ++l) { - const Eigen::Index sig = map3[l]; - ao(mu, nu, lam, sig) += cijk * c3[l] * prim(i, j, k, l); + const auto prim_data = buffer::get_raw_data(prim_T.buffer()); + + auto gamma_bra = detail_::gamma_ij(bs0, bs1); + auto gamma_ket = detail_::gamma_ij(bs2, bs3); + auto Q_bra = detail_::fine_k_ij(bs0, bs1); + auto Q_ket = detail_::fine_k_ij(bs2, bs3); + std::vector ao_data(naos[0] * naos[1] * naos[2] * naos[3], 0.0); + + auto ao_offset = [&](std::size_t i, std::size_t j, std::size_t k, + std::size_t l) { + return i * naos[1] * naos[2] * naos[3] + j * naos[2] * naos[3] + + k * naos[3] + l; + }; + auto prim_offset = [&](std::size_t i, std::size_t j, std::size_t k, + std::size_t l) { + return i * nprims[1] * nprims[2] * nprims[3] + + j * nprims[2] * nprims[3] + k * nprims[3] + l; + }; + + for(std::size_t pshell_i = 0; pshell_i < nprims[0]; ++pshell_i) { + const auto ci = c0[pshell_i]; + const auto mu = map0[pshell_i]; + const auto pi = pmap0[pshell_i]; + + for(std::size_t pshell_j = 0; pshell_j < nprims[1]; ++pshell_j) { + const auto nu = map1[pshell_j]; + const auto pj = pmap1[pshell_j]; + + const auto cij = ci * c1[pshell_j]; + const auto K_ij = K_bra[pi][pj]; + const auto gamma_ij = gamma_bra[pi][pj]; + const auto Q_ij = Q_bra[pi][pj]; + + for(std::size_t pshell_k = 0; pshell_k < nprims[2]; ++pshell_k) { + const auto lam = map2[pshell_k]; + const auto pk = pmap2[pshell_k]; + const double cijk = cij * c2[pshell_k]; + + for(std::size_t pshell_l = 0; pshell_l < nprims[3]; + ++pshell_l) { + const auto cl = c3[pshell_l]; + const auto sig = map3[pshell_l]; + const auto pl = pmap3[pshell_l]; + + const double K_kl = k_ket[pk][pl]; + const double Q_kl = Q_ket[pk][pl]; + const double gamma_kl = gamma_ket[pk][pl]; + + if(K_kl < thresh) continue; + if(K_ij * K_kl <= thresh) continue; + + auto Q_ijkl = Q_ij * Q_kl; + auto gamma_ijkl = gamma_ij + gamma_kl; + auto pfac = std::abs(Q_ijkl / std::sqrt(gamma_ijkl)); + if(pfac < thresh) continue; + + ao_data[ao_offset(mu, nu, lam, sig)] += + cijk * cl * + prim_data[prim_offset(pshell_i, pshell_j, pshell_k, + pshell_l)]; } } } } - // Wrap result in a TensorWrapper tensor - tensorwrapper::shape::Smooth shape( - {static_cast(n0), static_cast(n1), - static_cast(n2), static_cast(n3)}); + tensorwrapper::shape::Smooth shape({naos[0], naos[1], naos[2], naos[3]}); tensorwrapper::buffer::Contiguous buf(std::move(ao_data), shape); simde::type::tensor t(shape, std::move(buf)); diff --git a/src/integrals/utils/primitive_index_helpers.hpp b/src/integrals/utils/primitive_index_helpers.hpp new file mode 100644 index 00000000..3a75dbb9 --- /dev/null +++ b/src/integrals/utils/primitive_index_helpers.hpp @@ -0,0 +1,78 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include + +namespace integrals::utils { + +/** @brief Build a map from a primitive AO index to a primitive shell index. + * + * A contracted shell is a linear combination of primitive shells. Similarly, + * a contracted AO is a linear combination of primitive AOs. Given the offset + * of a primitive AO, this function returns the offset of the corresponding + * primitive shell. + * + * @param[in] bs The basis set to build the map for. + * @return A vector whose length is the number of primitive AO indices in the + * basis set. The value at index i is the primitive shell index that the i-th + * primitive AO belongs to. + */ +inline std::vector build_prim_ao_to_prim_shell_map( + const simde::type::ao_basis_set& bs) { + std::vector map; + std::size_t prim_offset = 0; + for(std::size_t s = 0; s < bs.n_shells(); ++s) { + const auto& shell = bs.shell(s); + const auto n_prims = shell.n_primitives(); + const auto n_aos = shell.size(); + for(std::size_t p = 0; p < n_prims; ++p) { + for(std::size_t a = 0; a < n_aos; ++a) { + map.push_back(prim_offset + p); + } + } + prim_offset += n_prims; + } + return map; +} + +/** @brief Maps a primitive AO index to a contracted AO index. + * + * @param[in] bs The basis set to build the map for. + * @return A vector whose length is the number of primitive AO indices in the + * basis set. The value at index i is the index of the contracted AO that + * the i-th primitive AO belongs to. + */ +inline std::vector build_prim_ao_to_cgto_map( + const simde::type::ao_basis_set& bs) { + std::vector map; + std::size_t contracted_ao = 0; + for(std::size_t s = 0; s < bs.n_shells(); ++s) { + const auto& shell = bs.shell(s); + const auto n_prims = shell.n_primitives(); + const auto n_aos = shell.size(); + for(std::size_t p = 0; p < n_prims; ++p) { + for(std::size_t a = 0; a < n_aos; ++a) { + map.push_back(contracted_ao + a); + } + } + contracted_ao += n_aos; + } + return map; +} + +} // namespace integrals::utils diff --git a/tests/cxx/unit/integrals/libint/detail_/test_primitive_pair_estimators.cpp b/tests/cxx/unit/integrals/libint/detail_/test_primitive_pair_estimators.cpp new file mode 100644 index 00000000..b903f0a7 --- /dev/null +++ b/tests/cxx/unit/integrals/libint/detail_/test_primitive_pair_estimators.cpp @@ -0,0 +1,86 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "integrals/libint/detail_/primitive_pair_estimators.hpp" +#undef DEPRECATED +#include "../../testing/testing.hpp" + +// Expected values below were recorded from a single run of each estimator on +// water_sto3g_basis_set() (same basis for both arguments). If the test basis or +// libint embedding changes, re-record by temporarily printing the outputs. + +using namespace integrals::testing; +using integrals::libint::detail_::coarse_k_ij; +using integrals::libint::detail_::fine_k_ij; +using integrals::libint::detail_::gamma_ij; +using integrals::libint::detail_::k_ij; + +namespace { + +constexpr double tol = 1.0e-11; + +} // namespace + +TEST_CASE("primitive_pair_estimators: gamma_ij") { + const auto aobs = water_sto3g_basis_set(); + REQUIRE(aobs.n_primitives() == 15); + const auto g = gamma_ij(aobs, aobs); + REQUIRE(g.size() == 15); + REQUIRE(g[0].size() == 15); + REQUIRE(g[0][0] == Catch::Approx(2.61418639999999982e+02).margin(tol)); + REQUIRE(g[0][14] == Catch::Approx(1.30878175403999990e+02).margin(tol)); + REQUIRE(g[14][0] == Catch::Approx(1.30878175403999990e+02).margin(tol)); + REQUIRE(g[14][14] == Catch::Approx(3.37710807999999973e-01).margin(tol)); + REQUIRE(g[1][2] == Catch::Approx(3.02524693000000013e+01).margin(tol)); +} + +TEST_CASE("primitive_pair_estimators: k_ij") { + const auto aobs = water_sto3g_basis_set(); + REQUIRE(aobs.n_primitives() == 15); + const auto K = k_ij(aobs, aobs); + REQUIRE(K.size() == 15); + REQUIRE(K[0].size() == 15); + for(std::size_t i = 0; i < 15; ++i) { + REQUIRE(K[i][i] == Catch::Approx(1.0).margin(1.0e-14)); + } + REQUIRE(K[0][9] == Catch::Approx(5.44974713611047489e-07).margin(1.0e-18)); +} + +TEST_CASE("primitive_pair_estimators: coarse_k_ij") { + const auto aobs = water_sto3g_basis_set(); + REQUIRE(aobs.n_primitives() == 15); + const auto Kc = coarse_k_ij(aobs, aobs); + REQUIRE(Kc.size() == 15); + REQUIRE(Kc[0].size() == 15); + REQUIRE(Kc[0][0] == Catch::Approx(1.80790216813702180e+01).margin(tol)); + REQUIRE(Kc[0][14] == Catch::Approx(1.71267463425960637e-01).margin(tol)); + REQUIRE(Kc[14][0] == Catch::Approx(1.71267463425960637e-01).margin(tol)); + REQUIRE(Kc[14][14] == Catch::Approx(6.96785377189288076e-03).margin(tol)); + REQUIRE(Kc[1][2] == Catch::Approx(5.27040829013411383e+00).margin(tol)); +} + +TEST_CASE("primitive_pair_estimators: fine_k_ij") { + const auto aobs = water_sto3g_basis_set(); + REQUIRE(aobs.n_primitives() == 15); + const auto Kf = fine_k_ij(aobs, aobs); + REQUIRE(Kf.size() == 15); + REQUIRE(Kf[0].size() == 15); + REQUIRE(Kf[0][0] == Catch::Approx(4.09063484384912246e-01).margin(tol)); + REQUIRE(Kf[0][14] == Catch::Approx(7.74033883652055620e-03).margin(tol)); + REQUIRE(Kf[14][0] == Catch::Approx(7.74033883652055620e-03).margin(tol)); + REQUIRE(Kf[14][14] == Catch::Approx(1.22041182423709954e-01).margin(tol)); + REQUIRE(Kf[1][2] == Catch::Approx(1.03047099111916607e+00).margin(tol)); +} diff --git a/tests/cxx/unit/integrals/libint/test_primitive_contractor.cpp b/tests/cxx/unit/integrals/libint/test_primitive_contractor.cpp index a0ff7686..4dcd0363 100644 --- a/tests/cxx/unit/integrals/libint/test_primitive_contractor.cpp +++ b/tests/cxx/unit/integrals/libint/test_primitive_contractor.cpp @@ -26,25 +26,45 @@ TEST_CASE("Primitive Contractor ERI4") { REQUIRE(mm.count("Primitive Contractor ERI4")); REQUIRE(mm.count("ERI4")); - auto aobs = water_sto3g_basis_set(); - for(const auto& type : - {chemist::ShellType::cartesian, chemist::ShellType::pure}) { - for(const auto& l : {1, 2, 3}) { - aobs.shell(2).l() = l; - aobs.shell(2).pure() = type; - simde::type::aos aos(aobs); - simde::type::aos_squared aos_squared(aos, aos); + auto corr_mod = mm.at("ERI4"); + auto test_mod = mm.at("Primitive Contractor ERI4"); - simde::type::v_ee_type op{}; + auto aobs = water_sto3g_basis_set(); + const auto is_cartesian = chemist::ShellType::cartesian; + const auto is_pure = chemist::ShellType::pure; + for(const auto& type : {is_cartesian, is_pure}) { + const auto type_str = + (type == is_cartesian) ? "Cartesian" : "Spherical"; + SECTION(type_str) { + for(const auto& l : {1, 2, 3}) { + SECTION("Angular Momentum: " + std::to_string(l)) { + aobs.shell(2).l() = l; + aobs.shell(2).pure() = type; + simde::type::aos aos(aobs); + simde::type::aos_squared aos_squared(aos, aos); + simde::type::v_ee_type op{}; + chemist::braket::BraKet braket(aos_squared, op, + aos_squared); - chemist::braket::BraKet braket(aos_squared, op, aos_squared); - auto corr_mod = mm.at("ERI4"); - auto test_mod = mm.at("Primitive Contractor ERI4"); + for(const auto& thresh : {1e-16, 1e-8, 1e-6}) { + SECTION("Screening Threshold: " + + std::to_string(thresh)) { + auto test_mod_copy = test_mod.unlocked_copy(); + test_mod_copy.change_input("Screening Threshold", + thresh); - auto T_contracted = test_mod.run_as(braket); - auto T_eri4 = corr_mod.run_as(braket); + auto corr_mod_copy = corr_mod.unlocked_copy(); + corr_mod_copy.change_input("Threshold", thresh); + auto T_contracted = + test_mod_copy.run_as(braket); + auto T_eri4 = corr_mod_copy.run_as(braket); - REQUIRE(approximately_equal(T_contracted, T_eri4, 1.0e-14)); + REQUIRE(approximately_equal(T_contracted, T_eri4, + 1.0e-14)); + } + } + } + } } } }