From 7a5d9eb250f2b0e9e96be1bb672776f15a833565 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 6 Feb 2026 15:18:48 -0600 Subject: [PATCH 01/14] backup --- include/integrals/integrals.hpp | 3 +- include/integrals/integrals_mm.hpp | 2 + include/integrals/property_types.hpp | 40 ++++ src/integrals/ao_integrals/ao_integrals.hpp | 1 + src/integrals/ao_integrals/uq_driver.cpp | 46 ++-- src/integrals/integrals_mm.cpp | 4 +- src/integrals/libint/analytic_error.cpp | 40 ++++ .../libint/black_box_primitive_estimator.cpp | 78 +++++++ .../cauchy_schwarz_primitive_estimator.cpp | 114 ++++++++++ src/integrals/libint/libint.cpp | 13 +- src/integrals/libint/libint.hpp | 19 +- .../libint/primitive_error_model.cpp | 202 ++++++++++++++++++ src/integrals/utils/rank2_shell_norm.hpp | 41 ++++ .../ao_integrals/ao_integrals_driver.cpp | 1 + .../integrals/ao_integrals/coulomb_metric.cpp | 1 + .../integrals/ao_integrals/df_integral.cpp | 1 + .../ao_integrals/j_density_fitted.cpp | 1 + .../integrals/ao_integrals/j_four_center.cpp | 1 + .../ao_integrals/k_density_fitted.cpp | 1 + .../integrals/ao_integrals/k_four_center.cpp | 1 + .../integrals/ao_integrals/test_uq_driver.cpp | 3 +- .../unit/integrals/libint/analytic_error.cpp | 34 +++ .../libint/black_box_primitive_estimator.cpp | 62 ++++++ .../cauchy_schwarz_primitive_estimator.cpp | 31 +++ .../libint/primitive_error_model.cpp | 41 ++++ .../cxx/unit/integrals/libint/test_error.hpp | 77 +++++++ tests/cxx/unit/integrals/testing.hpp | 7 + .../unit/integrals/utils/rank2_shell_norm.cpp | 99 +++++++++ 28 files changed, 924 insertions(+), 40 deletions(-) create mode 100644 src/integrals/libint/analytic_error.cpp create mode 100644 src/integrals/libint/black_box_primitive_estimator.cpp create mode 100644 src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp create mode 100644 src/integrals/libint/primitive_error_model.cpp create mode 100644 src/integrals/utils/rank2_shell_norm.hpp create mode 100644 tests/cxx/unit/integrals/libint/analytic_error.cpp create mode 100644 tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp create mode 100644 tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp create mode 100644 tests/cxx/unit/integrals/libint/primitive_error_model.cpp create mode 100644 tests/cxx/unit/integrals/libint/test_error.hpp create mode 100644 tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp diff --git a/include/integrals/integrals.hpp b/include/integrals/integrals.hpp index ac432349..e1d2e352 100644 --- a/include/integrals/integrals.hpp +++ b/include/integrals/integrals.hpp @@ -22,7 +22,8 @@ * your unit test also needs most of the headers included by it). */ #pragma once -#include "integrals/integrals_mm.hpp" +#include +#include /** @namespace integrals * diff --git a/include/integrals/integrals_mm.hpp b/include/integrals/integrals_mm.hpp index 55e479d7..28c6ae9a 100644 --- a/include/integrals/integrals_mm.hpp +++ b/include/integrals/integrals_mm.hpp @@ -27,4 +27,6 @@ namespace integrals { */ DECLARE_PLUGIN(integrals); +void set_defaults(pluginplay::ModuleManager&); + } // end namespace integrals diff --git a/include/integrals/property_types.hpp b/include/integrals/property_types.hpp index aa275d5a..2798c44f 100644 --- a/include/integrals/property_types.hpp +++ b/include/integrals/property_types.hpp @@ -30,6 +30,46 @@ */ namespace integrals::property_types { +// PT used to estimate the contribution of primitive pairs +DECLARE_PROPERTY_TYPE(PrimitivePairEstimator); +PROPERTY_TYPE_INPUTS(PrimitivePairEstimator) { + using ao_basis = const simde::type::ao_basis_set&; + auto rv = pluginplay::declare_input() + .add_field("Bra Basis Set") + .add_field("Ket Basis Set"); + rv["Bra Basis Set"].set_description( + "The atomic orbital basis set for the bra"); + rv["Ket Basis Set"].set_description( + "The atomic orbital basis set for the ket"); + return rv; +} + +PROPERTY_TYPE_RESULTS(PrimitivePairEstimator) { + using tensor = simde::type::tensor; + auto rv = pluginplay::declare_result().add_field( + "Primitive Pair Estimates"); + rv["Primitive Pair Estimates"].set_description( + "A tensor containing the estimated values for each primitive pair " + "integral"); + return rv; +} + +template +DECLARE_TEMPLATED_PROPERTY_TYPE(Uncertainty, BasePT); + +template +TEMPLATED_PROPERTY_TYPE_INPUTS(Uncertainty, BasePT) { + auto rv = BasePT::inputs(); + auto rv0 = rv.template add_field("Tolerance"); + rv0["Tolerance"].set_description("The screening threshold"); + return rv0; +} + +template +TEMPLATED_PROPERTY_TYPE_RESULTS(Uncertainty, BasePT) { + return BasePT::results(); +} + using DecontractBasisSet = simde::Convert; diff --git a/src/integrals/ao_integrals/ao_integrals.hpp b/src/integrals/ao_integrals/ao_integrals.hpp index de519eed..cee09ef4 100644 --- a/src/integrals/ao_integrals/ao_integrals.hpp +++ b/src/integrals/ao_integrals/ao_integrals.hpp @@ -40,6 +40,7 @@ inline void set_defaults(pluginplay::ModuleManager& mm) { mm.change_submod("Density Fitting Integral", "Coulomb Metric", "Coulomb Metric"); mm.change_submod("UQ Driver", "ERIs", "ERI4"); + mm.change_submod("UQ Driver", "ERI Error", "Analytic Error"); } inline void load_modules(pluginplay::ModuleManager& mm) { diff --git a/src/integrals/ao_integrals/uq_driver.cpp b/src/integrals/ao_integrals/uq_driver.cpp index 5b309a85..9d03ba0c 100644 --- a/src/integrals/ao_integrals/uq_driver.cpp +++ b/src/integrals/ao_integrals/uq_driver.cpp @@ -15,6 +15,7 @@ */ #include "ao_integrals.hpp" +#include using namespace tensorwrapper; @@ -37,19 +38,21 @@ struct Kernel { const std::span error) { Tensor rv; - if constexpr(types::is_uncertain_v) { - auto rv_buffer = buffer::make_contiguous(m_shape); - auto rv_data = buffer::get_raw_data(rv_buffer); + using float_type = std::decay_t; + if constexpr(types::is_uncertain_v) { + throw std::runtime_error("Did not expect an uncertain type"); + } else { + using uq_type = sigma::Uncertain; + auto rv_buffer = buffer::make_contiguous(m_shape); + auto rv_data = buffer::get_raw_data(rv_buffer); for(std::size_t i = 0; i < t.size(); ++i) { - const auto elem = t[i].mean(); - const auto elem_error = error[i].mean(); - rv_data[i] = FloatType(elem, elem_error); + const auto elem = t[i]; + const auto elem_error = error[i]; + rv_data[i] = uq_type(elem, elem_error); } - rv = tensorwrapper::Tensor(m_shape, std::move(rv_buffer)); - } else { - throw std::runtime_error("Expected an uncertain type"); } + return rv; } shape_type m_shape; @@ -63,35 +66,24 @@ UQ Integrals Driver } // namespace -using eri_pt = simde::ERI4; +using eri_pt = simde::ERI4; +using error_pt = integrals::property_types::Uncertainty; MODULE_CTOR(UQDriver) { satisfies_property_type(); description(desc); add_submodule("ERIs"); - add_input("benchmark precision").set_default(1.0e-16); - add_input("precision").set_default(1.0e-16); + add_submodule("ERI Error"); } MODULE_RUN(UQDriver) { - auto tau_0 = inputs.at("benchmark precision").value(); - auto tau = inputs.at("precision").value(); + const auto& [braket] = eri_pt::unwrap_inputs(inputs); auto& eri_mod = submods.at("ERIs").value(); + auto tol = eri_mod.inputs().at("Threshold").value(); - auto benchmark_mod = eri_mod.unlocked_copy(); - benchmark_mod.change_input("Threshold", tau_0); - benchmark_mod.change_input("With UQ?", true); - - auto normal_mod = eri_mod.unlocked_copy(); - normal_mod.change_input("Threshold", tau); - normal_mod.change_input("With UQ?", true); - - const auto& [t_0] = eri_pt::unwrap_results(benchmark_mod.run(inputs)); - const auto& [t] = eri_pt::unwrap_results(normal_mod.run(inputs)); - - simde::type::tensor error; - error("m,n,l,s") = t("m,n,l,s") - t_0("m,n,l,s"); + const auto& t = eri_mod.run_as(braket); + const auto& error = submods.at("ERI Error").run_as(braket, tol); using buffer::visit_contiguous_buffer; shape::Smooth shape = t.buffer().layout().shape().as_smooth().make_smooth(); diff --git a/src/integrals/integrals_mm.cpp b/src/integrals/integrals_mm.cpp index 414b7923..7021528c 100644 --- a/src/integrals/integrals_mm.cpp +++ b/src/integrals/integrals_mm.cpp @@ -28,6 +28,8 @@ namespace integrals { * @throw none No throw guarantee */ void set_defaults(pluginplay::ModuleManager& mm) { + libint::set_defaults(mm); + ao_integrals::set_defaults(mm); mm.change_submod("AO integral driver", "Kinetic", "Kinetic"); mm.change_submod("AO integral driver", "Electron-Nuclear attraction", "Nuclear"); @@ -41,8 +43,6 @@ void load_modules(pluginplay::ModuleManager& mm) { ao_integrals::load_modules(mm); libint::load_modules(mm); utils::load_modules(mm); - set_defaults(mm); - ao_integrals::set_defaults(mm); } } // namespace integrals diff --git a/src/integrals/libint/analytic_error.cpp b/src/integrals/libint/analytic_error.cpp new file mode 100644 index 00000000..e4b01005 --- /dev/null +++ b/src/integrals/libint/analytic_error.cpp @@ -0,0 +1,40 @@ +#include "libint.hpp" +#include + +namespace integrals::libint { +namespace { + +const auto desc = "Uses the error in the ERI4 as the uncertainty."; +} + +using eri4_pt = simde::ERI4; +using pt = integrals::property_types::Uncertainty; + +MODULE_CTOR(AnalyticError) { + satisfies_property_type(); + description(desc); + + add_submodule("ERI4s"); +} + +MODULE_RUN(AnalyticError) { + const auto& [braket, tol] = pt::unwrap_inputs(inputs); + + auto& eri_mod = submods.at("ERI4s"); + + auto normal_mod = eri_mod.value().unlocked_copy(); + normal_mod.change_input("Threshold", tol); + + // N.b., t_0 is the benchmark value + const auto& t_0 = eri_mod.run_as(braket); + const auto& t = normal_mod.run_as(braket); + + simde::type::tensor error; + error("m,n,l,s") = t("m,n,l,s") - t_0("m,n,l,s"); + + // Wrap and return the results + auto rv = results(); + return pt::wrap_results(rv, error); +} + +} // namespace integrals::libint \ No newline at end of file diff --git a/src/integrals/libint/black_box_primitive_estimator.cpp b/src/integrals/libint/black_box_primitive_estimator.cpp new file mode 100644 index 00000000..4a72119c --- /dev/null +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -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. + */ + +#include "libint.hpp" +#include +#include + +namespace integrals::libint { +namespace { + +const auto desc = ""; + +} + +using pt = integrals::property_types::PrimitivePairEstimator; + +MODULE_CTOR(BlackBoxPrimitiveEstimator) { + satisfies_property_type(); + description(desc); + // TODO: Add citation for Chemist paper +} + +MODULE_RUN(BlackBoxPrimitiveEstimator) { + const auto&& [bra_basis, ket_basis] = pt::unwrap_inputs(inputs); + + const auto n_bra = bra_basis.n_primitives(); + const auto n_ket = ket_basis.n_primitives(); + + // TODO: Use floating point type of the basis sets + using float_type = double; + std::vector buffer(n_bra * n_ket, 0.0); + + using iter_type = std::decay_t; // Type of the loop index + + for(iter_type bra_i = 0; bra_i < n_bra; ++bra_i) { + const auto bra_prim_i = bra_basis.primitive(bra_i); + const auto bra_zeta = bra_prim_i.exponent(); + const auto bra_coeff = std::fabs(bra_prim_i.coefficient()); + const auto bra_center = bra_prim_i.center(); + const auto bra_offset = bra_i * n_ket; + + for(iter_type ket_i = 0; ket_i < n_ket; ++ket_i) { + const auto ket_prim_i = ket_basis.primitive(ket_i); + const auto ket_zeta = ket_prim_i.exponent(); + const auto ket_coeff = std::fabs(ket_prim_i.coefficient()); + const auto ket_center = ket_prim_i.center(); + + // This is "K bar" in Eq. 11 in the SI of the Chemist paper + const auto dist = (bra_center - ket_center).magnitude(); + const auto dist2 = dist * dist; + const auto ratio = bra_zeta * ket_zeta / (bra_zeta + ket_zeta); + const auto coeff = bra_coeff * ket_coeff; + buffer[bra_offset + ket_i] = coeff * std::exp(-1.0 * ratio * dist2); + } + } + + tensorwrapper::shape::Smooth shape({n_bra, n_ket}); + tensorwrapper::buffer::Contiguous tw_buffer(std::move(buffer), shape); + simde::type::tensor rv(shape, std::move(tw_buffer)); + + auto result = results(); + return pt::wrap_results(result, rv); +} + +} // namespace integrals::libint \ No newline at end of file diff --git a/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp new file mode 100644 index 00000000..a312d947 --- /dev/null +++ b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -0,0 +1,114 @@ +/* + * 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 "../utils/rank2_shell_norm.hpp" +#include "libint.hpp" +#include +#include + +namespace integrals::libint { +namespace { + +const auto desc = ""; + +} // namespace + +using decontract_pt = integrals::property_types::DecontractBasisSet; +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"); +} + +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 = bra_prims.n_primitives(); + const auto n_ket = ket_prims.n_primitives(); + + // Should always be true, but we check for sanity + assert(n_bra == bra_basis.n_primitives()); + assert(n_ket == ket_basis.n_primitives()); + + // TODO: Only get the hyper diagonal we actually need + simde::type::aos_squared bra(bra_prims, ket_prims); + simde::type::v_ee_type v_ee{}; + chemist::braket::BraKet mnls(bra, v_ee, bra); + const auto& prim4 = submods.at("ERI4").run_as(mnls); + + 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 buffer(n_bra * n_ket, 0.0); + + using iter_type = std::decay_t; // Type of the loop index + + // We loop over the decontracted basis, but make sure to grab the + // contraction coefficients from the real basis set. By looping over the + // decontracted basis we can make this look like the "usual" shell pair + // loop + // N.b. there's one shell per primitive in the real basis + // N.b. since the basis is decontracted the number of Cartesian/spherical + // AOs in a shell is the same as the number of Cartesian/spherical + // primitives in that shell + + std::array offsets{0, 0, 0, 0}; + std::array naos{0, 0, 0, 0}; + for(iter_type shell_i = 0; shell_i < bra_prims.n_shells(); ++shell_i) { + const auto bra_prim_i = bra_basis.primitive(shell_i); + const auto bra_coeff = std::fabs(bra_prim_i.coefficient()); + naos[0] = bra_prims.shell(shell_i).size(); + naos[2] = naos[0]; + + const auto bra_shell_offset = shell_i * n_ket; + offsets[1] = 0; + offsets[3] = 0; + for(iter_type shell_j = 0; shell_j < ket_prims.n_shells(); ++shell_j) { + const auto ket_prim_i = ket_basis.primitive(shell_j); + const auto ket_coeff = std::fabs(ket_prim_i.coefficient()); + naos[1] = ket_prims.shell(shell_j).size(); + naos[3] = naos[1]; + + auto C_ab = bra_coeff * ket_coeff; + auto shell_norm = utils::rank2_shell_norm(eris, offsets, naos); + buffer[bra_shell_offset + shell_j] = C_ab * std::sqrt(shell_norm); + + offsets[1] += naos[1]; + offsets[3] += naos[3]; + } + offsets[0] += naos[0]; + offsets[2] += naos[2]; + } + + tensorwrapper::shape::Smooth shape({n_bra, n_ket}); + tensorwrapper::buffer::Contiguous tw_buffer(std::move(buffer), shape); + simde::type::tensor rv(shape, std::move(tw_buffer)); + + auto result = results(); + return pt::wrap_results(result, rv); +} + +} // namespace integrals::libint \ No newline at end of file diff --git a/src/integrals/libint/libint.cpp b/src/integrals/libint/libint.cpp index d391c2a0..4395a1e5 100644 --- a/src/integrals/libint/libint.cpp +++ b/src/integrals/libint/libint.cpp @@ -192,7 +192,12 @@ EXTERN_LIBINT(aos_squared, v_ee_type, aos_squared); #undef EXTERN_LIBINT void set_defaults(pluginplay::ModuleManager& mm) { - // Set any default associations + mm.change_submod("Primitive Error Model", "Primitive Pair Estimator", + "Black Box Primitive Pair Estimator"); + mm.change_submod("CauchySchwarz Estimator", "Decontract Basis Set", + "Decontract Basis Set"); + mm.change_submod("CauchySchwarz Estimator", "ERI4", "ERI4"); + mm.change_submod("Analytic Error", "ERI4s", "ERI4"); } #define LOAD_LIBINT(bra, op, ket, key) mm.add_module(key) @@ -208,7 +213,11 @@ void load_modules(pluginplay::ModuleManager& mm) { LOAD_LIBINT(aos, v_ee_type, aos, "ERI2"); LOAD_LIBINT(aos, v_ee_type, aos_squared, "ERI3"); LOAD_LIBINT(aos_squared, v_ee_type, aos_squared, "ERI4"); - set_defaults(mm); + mm.add_module( + "Black Box Primitive Pair Estimator"); + mm.add_module("CauchySchwarz Estimator"); + mm.add_module("Primitive Error Model"); + mm.add_module("Analytic Error"); } #undef LOAD_LIBINT diff --git a/src/integrals/libint/libint.hpp b/src/integrals/libint/libint.hpp index d05968fb..46fcfa04 100644 --- a/src/integrals/libint/libint.hpp +++ b/src/integrals/libint/libint.hpp @@ -24,6 +24,18 @@ */ namespace integrals::libint { +/** @brief The Module for computing AO Integrals + * + * @tparam BraKetType The type of the BraKet input + */ +template +DECLARE_MODULE(Libint); + +DECLARE_MODULE(BlackBoxPrimitiveEstimator); +DECLARE_MODULE(CauchySchwarzPrimitiveEstimator); +DECLARE_MODULE(PrimitiveErrorModel); +DECLARE_MODULE(AnalyticError); + using simde::type::braket; using simde::type::aos; @@ -35,13 +47,6 @@ using simde::type::t_e_type; using simde::type::v_ee_type; using simde::type::v_en_type; -/** @brief The Module for computing AO Integrals - * - * @tparam BraKetType The type of the BraKet input - */ -template -DECLARE_MODULE(Libint); - /** @brief Load the libint modules into a Module Manager * * @param mm The Module Manager to load the modules into diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp new file mode 100644 index 00000000..5ac49c0a --- /dev/null +++ b/src/integrals/libint/primitive_error_model.cpp @@ -0,0 +1,202 @@ +#include "libint.hpp" +#include +#include + +namespace integrals::libint { +namespace { + +const auto desc = "Compute uncertainty estimates for ERI4 integrals"; + +// Sums contributions Q_abQ_cd when Q_ab or Q_cd is below tol +template +auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, +TensorType&& Q_ab, TensorType&& Q_cd, double tol){ + // Get the number of primitives in each shell + const auto n_prim0 = shells[0].n_primitives(); + const auto n_prim1 = shells[1].n_primitives(); + const auto n_prim2 = shells[2].n_primitives(); + const auto n_prim3 = shells[3].n_primitives(); + + using float_type = double; // TODO: get from Q_ab or Q_cd + using iter_type = std::decay_t; // Type of the loop index + + // Create an array to hold the primitive indices + std::array prim{0, 0, 0, 0}; + + // This is the accumulated error for this shell block + double error = 0.0; + + + using wtf::fp::float_cast; + for(prim[0] = 0; prim[0] < n_prim0; ++prim[0]){ + const auto p0 = prim_offsets[0] + prim[0]; + + for(prim[1] = 0; prim[1] < n_prim1; ++prim[1]){ + const auto p1 = prim_offsets[1] + prim[1]; + + // Was the Q_ab element neglected for primitive pair (p0, p1)? + std::vector i01{p0, p1}; + + auto Q_01 = float_cast(Q_ab.get_elem(i01)); + auto Q_01_mag = std::fabs(Q_01); + + for(prim[2] = 0; prim[2] < n_prim2; ++prim[2]){ + const auto p2 = prim_offsets[2] + prim[2]; + + for(prim[3] = 0; prim[3] < n_prim3; ++prim[3]){ + const auto p3 = prim_offsets[3] + prim[3]; + + // Was the Q_cd element neglected for pair (p2, p3)? + std::vector i23{p2, p3}; + auto Q_23 = float_cast(Q_cd.get_elem(i23)); + auto Q_23_mag = std::fabs(Q_23); + + // If either was neglected we pick up an error Q_01 * Q_23 + if(Q_01_mag * Q_23_mag < tol) + error += Q_01_mag * Q_23_mag; + } // End loop over prim[3] + } // End loop over prim[2] + } // End loop over prim[1] + } // End loop over prim[0] + + return error; +} + +// Sets all AO integrals in the shell block to the computed error +template +void fill_ao_block(ShellsType&& shells, OffsetTypes&& ao_offsets, TensorType&& t, + double error){ + + // Get the number of AOs in each shell + const auto n_aos0 = shells[0].size(); + const auto n_aos1 = shells[1].size(); + const auto n_aos2 = shells[2].size(); + const auto n_aos3 = shells[3].size(); + + // Make an array to hold the AO indices + using iter_type = std::decay_t; // Type of the loop index + std::array ao{0, 0, 0, 0}; + + for(ao[0] = 0; ao[0] < n_aos0; ++ao[0]){ + const auto a0 = ao_offsets[0] + ao[0]; + + for(ao[1] = 0; ao[1] < n_aos1; ++ao[1]){ + const auto a1 = ao_offsets[1] + ao[1]; + + for(ao[2] = 0; ao[2] < n_aos2; ++ao[2]){ + const auto a2 = ao_offsets[2] + ao[2]; + + for(ao[3] = 0; ao[3] < n_aos3; ++ao[3]){ + const auto a3 = ao_offsets[3] + ao[3]; + + std::vector index{a0, a1, a2, a3}; + t.set_elem(index, error); + } // End loop over ao[3] + } // End loop over ao[2] + } // End loop over ao[1] + } // End loop over ao[0] +} + +} // namespace + +using eri4_pt = simde::ERI4; +using pair_estimator_pt = integrals::property_types::PrimitivePairEstimator; +using pt = integrals::property_types::Uncertainty; + +MODULE_CTOR(PrimitiveErrorModel) { + satisfies_property_type(); + description(desc); + // TODO citation for Chemist paper + + add_submodule("Primitive Pair Estimator") + .set_description("The module used to estimate the contributions of " + "primitive pairs to the overall integral values"); +} + +MODULE_RUN(PrimitiveErrorModel) { + const auto& [braket, tol] = pt::unwrap_inputs(inputs); + + const auto bra0 = braket.bra().first.ao_basis_set(); + const auto bra1 = braket.bra().second.ao_basis_set(); + const auto ket0 = braket.ket().first.ao_basis_set(); + const auto ket1 = braket.ket().second.ao_basis_set(); + + // Get the Q_ab and Q_cd tensors + auto& estimator = submods.at("Primitive Pair Estimator"); + const auto& Q_ab_tw = estimator.run_as(bra0, bra1); + const auto& Q_cd_tw = estimator.run_as(ket0, ket1); + + // XXX: Workaround for needing the contiguous buffers to access elements + using tensorwrapper::buffer::make_contiguous; + const auto& Q_ab = make_contiguous(Q_ab_tw.buffer()); + const auto& Q_cd = make_contiguous(Q_cd_tw.buffer()); + + // Get the number of AOs in each basis set + const auto n_mu = bra0.n_aos(); + const auto n_nu = bra1.n_aos(); + const auto n_lambda = ket0.n_aos(); + const auto n_sigma = ket1.n_aos(); + + // Make the return buffer + using float_type = double; // TODO: get from Q_ab or Q_cd + tensorwrapper::shape::Smooth shape({n_mu, n_nu, n_lambda, n_sigma}); + std::vector raw_buffer(shape.size(), 0.0); + tensorwrapper::buffer::Contiguous buffer(std::move(raw_buffer), shape); + + // Make arrays to hold the indices and offsets + using iter_type = std::decay_t; // Type of the loop index + using shell_view = decltype(bra0.shell(0)); + std::array n_shells{ + bra0.n_shells(), bra1.n_shells(), ket0.n_shells(), ket1.n_shells()}; + std::array shell_i{0, 0, 0, 0}; + std::array prim_offset{0, 0, 0, 0}; + std::array ao_offsets{0, 0, 0, 0}; + + // We'll collect the shells into this array + std::array shells; + + for(shell_i[0] = 0 ; shell_i[0] < n_shells[0]; ++shell_i[0]){ + shells[0] = bra0.shell(shell_i[0]); + + prim_offset[1] = 0; + ao_offsets[1] = 0; + for(shell_i[1] = 0; shell_i[1] < n_shells[1]; ++shell_i[1]){ + shells[1] = bra1.shell(shell_i[1]); + + prim_offset[2] = 0; + ao_offsets[2] = 0; + for(shell_i[2] = 0; shell_i[2] < n_shells[2]; ++shell_i[2]){ + shells[2] = ket0.shell(shell_i[2]); + + prim_offset[3] = 0; + ao_offsets[3] = 0; + for(shell_i[3] = 0; shell_i[3] < n_shells[3]; ++shell_i[3]){ + shells[3] = ket1.shell(shell_i[3]); + + auto shell_error = + shell_block_error(shells, prim_offset, Q_ab, Q_cd, tol); + fill_ao_block(shells, ao_offsets, buffer, shell_error); + + prim_offset[3] += shells[3].n_primitives(); + ao_offsets[3] += shells[3].size(); + } // End loop over shell_i[3] + + prim_offset[2] += shells[2].n_primitives(); + ao_offsets[2] += shells[2].size(); + }// End loop over shell_i[2] + + prim_offset[1] += shells[1].n_primitives(); + ao_offsets[1] += shells[1].size(); + } // End loop over shell_i[1] + + prim_offset[0] += shells[0].n_primitives(); + ao_offsets[0] += shells[0].size(); + } // End loop over shell_i[0] + + + simde::type::tensor error(shape, std::move(buffer)); + auto result = results(); + return pt::wrap_results(result, error); +} + +} \ No newline at end of file diff --git a/src/integrals/utils/rank2_shell_norm.hpp b/src/integrals/utils/rank2_shell_norm.hpp new file mode 100644 index 00000000..ab522569 --- /dev/null +++ b/src/integrals/utils/rank2_shell_norm.hpp @@ -0,0 +1,41 @@ +#pragma once +#include +#include +#include +#include + +namespace integrals::utils { + +template +auto rank2_shell_norm(TensorType&& buffer, std::array offsets, + std::array naos) { + using float_type = double; + // Take the shell norm over the primitive integrals in this shell pair. + float_type shell_norm = 0.0; + std::array ao{0, 0, 0, 0}; + std::vector abs_ao(4, 0); // Absolute indices + + using wtf::fp::float_cast; + for(ao[0] = 0; ao[0] < naos[0]; ++ao[0]) { + abs_ao[0] = offsets[0] + ao[0]; // Absolute index + + for(ao[1] = 0; ao[1] < naos[1]; ++ao[1]) { + abs_ao[1] = offsets[1] + ao[1]; + + for(ao[2] = 0; ao[2] < naos[2]; ++ao[2]) { + abs_ao[2] = offsets[2] + ao[2]; + + for(ao[3] = 0; ao[3] < naos[3]; ++ao[3]) { + abs_ao[3] = offsets[3] + ao[3]; + + const auto I_ijkl = buffer.get_elem(abs_ao); + auto as_float = std::fabs(float_cast(I_ijkl)); + shell_norm += as_float * as_float; + } // End loop over ao[3] + } // End loop over ao[2] + } // End loop over ao[1] + } // End loop over ao[0] + + return std::sqrt(shell_norm); +} +} // namespace integrals::utils \ No newline at end of file diff --git a/tests/cxx/unit/integrals/ao_integrals/ao_integrals_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/ao_integrals_driver.cpp index 479375d3..d3ac55f2 100644 --- a/tests/cxx/unit/integrals/ao_integrals/ao_integrals_driver.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/ao_integrals_driver.cpp @@ -53,6 +53,7 @@ TEST_CASE("AOIntegralsDriver") { pluginplay::ModuleManager mm; integrals::load_modules(mm); + integrals::set_defaults(mm); REQUIRE(mm.count("AO integral driver")); auto& mod = mm.at("AO integral driver"); diff --git a/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp b/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp index 6ceea87d..b49e4c48 100644 --- a/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp @@ -21,6 +21,7 @@ TEST_CASE("Coulomb Metric") { pluginplay::ModuleManager mm; integrals::load_modules(mm); + integrals::set_defaults(mm); REQUIRE(mm.count("Coulomb Metric")); // Get basis set diff --git a/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp b/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp index e3045bf4..f58bbd21 100644 --- a/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp @@ -21,6 +21,7 @@ TEST_CASE("Density Fitting Integral") { pluginplay::ModuleManager mm; integrals::load_modules(mm); + integrals::set_defaults(mm); REQUIRE(mm.count("Density Fitting Integral")); // Get basis set diff --git a/tests/cxx/unit/integrals/ao_integrals/j_density_fitted.cpp b/tests/cxx/unit/integrals/ao_integrals/j_density_fitted.cpp index 95777064..17938991 100644 --- a/tests/cxx/unit/integrals/ao_integrals/j_density_fitted.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/j_density_fitted.cpp @@ -21,6 +21,7 @@ TEST_CASE("Density Fitted J builder") { pluginplay::ModuleManager mm; integrals::load_modules(mm); + integrals::set_defaults(mm); REQUIRE(mm.count("Density Fitted J builder")); // Get basis set diff --git a/tests/cxx/unit/integrals/ao_integrals/j_four_center.cpp b/tests/cxx/unit/integrals/ao_integrals/j_four_center.cpp index c455391b..fe882f9b 100644 --- a/tests/cxx/unit/integrals/ao_integrals/j_four_center.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/j_four_center.cpp @@ -21,6 +21,7 @@ TEST_CASE("Four center J builder") { pluginplay::ModuleManager mm; integrals::load_modules(mm); + integrals::set_defaults(mm); REQUIRE(mm.count("Four center J builder")); // Get basis set diff --git a/tests/cxx/unit/integrals/ao_integrals/k_density_fitted.cpp b/tests/cxx/unit/integrals/ao_integrals/k_density_fitted.cpp index f2341506..b08664c7 100644 --- a/tests/cxx/unit/integrals/ao_integrals/k_density_fitted.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/k_density_fitted.cpp @@ -21,6 +21,7 @@ TEST_CASE("Density Fitted K builder") { pluginplay::ModuleManager mm; integrals::load_modules(mm); + integrals::set_defaults(mm); REQUIRE(mm.count("Density Fitted K builder")); // Get basis set diff --git a/tests/cxx/unit/integrals/ao_integrals/k_four_center.cpp b/tests/cxx/unit/integrals/ao_integrals/k_four_center.cpp index 97da531e..d5010ab2 100644 --- a/tests/cxx/unit/integrals/ao_integrals/k_four_center.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/k_four_center.cpp @@ -21,6 +21,7 @@ TEST_CASE("Four center K builder") { pluginplay::ModuleManager mm; integrals::load_modules(mm); + integrals::set_defaults(mm); REQUIRE(mm.count("Four center K builder")); // Get basis set diff --git a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp index 055f1280..e7a9b722 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp @@ -53,9 +53,10 @@ TEST_CASE("UQ Driver") { auto rt = std::make_unique(); pluginplay::ModuleManager mm(std::move(rt), nullptr); integrals::load_modules(mm); + integrals::set_defaults(mm); REQUIRE(mm.count("UQ Driver")); - mm.change_input("UQ Driver", "precision", 1.0e-6); + mm.change_input("ERI4", "Threshold", 1.0e-6); // Get basis set auto mol = test::h2_molecule(); diff --git a/tests/cxx/unit/integrals/libint/analytic_error.cpp b/tests/cxx/unit/integrals/libint/analytic_error.cpp new file mode 100644 index 00000000..1c95825f --- /dev/null +++ b/tests/cxx/unit/integrals/libint/analytic_error.cpp @@ -0,0 +1,34 @@ +#include "test_error.hpp" +#include + +using eri4_pt = simde::ERI4; +using pt = integrals::property_types::Uncertainty; + +using namespace integrals::libint::test; + +TEST_CASE("AnalyticError") { + pluginplay::ModuleManager mm; + integrals::load_modules(mm); + integrals::set_defaults(mm); + + auto& mod = mm.at("Analytic Error"); + + simde::type::v_ee_type v_ee{}; + using float_type = double; + + SECTION("H2 Dimer/STO-3G (03|12)") { + auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + + simde::type::aos_squared bra(bra0, bra1); + simde::type::aos_squared ket(ket0, ket1); + chemist::braket::BraKet mnls(bra, v_ee, ket); + float_type tol = 1.0e-10; + auto error = mod.run_as(mnls, tol); + tensorwrapper::shape::Smooth shape({1, 1, 1, 1}); + auto benchmark = test::make_tw_buffer(0.0002263495626484, shape); + auto screened = test::make_tw_buffer(0.0002263440759384, shape); + tensorwrapper::buffer::Contiguous corr(benchmark); + corr("m,n,l,s") = screened("m,n,l,s") - benchmark("m,n,l,s"); + REQUIRE(error.buffer().approximately_equal(corr, 1.0e-12)); + } +} \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp new file mode 100644 index 00000000..de2a07f6 --- /dev/null +++ b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp @@ -0,0 +1,62 @@ +/* + * 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 "testing.hpp" +#include + +using pt = integrals::property_types::PrimitivePairEstimator; + +namespace { + +std::vector h2_h2_corr{ + 0.023817430147884473, 0.08261663936778646, 0.06861998972363427, + 1.4555047643750448e-05, 0.00844595180076824, 0.03423471319097009, + 0.08261663936778646, 0.28657621993836907, 0.23802538347833696, + 0.00844595180076824, 0.07444365105885908, 0.13404298325096972, + 0.06861998972363427, 0.23802538347833696, 0.19769987611740356, + 0.03423471319097009, 0.13404298325096972, 0.1372685140907819, + 1.4555047643750448e-05, 0.00844595180076824, 0.03423471319097009, + 0.023817430147884473, 0.08261663936778646, 0.06861998972363427, + 0.00844595180076824, 0.07444365105885908, 0.13404298325096972, + 0.08261663936778646, 0.28657621993836907, 0.23802538347833696, + 0.03423471319097009, 0.13404298325096972, 0.1372685140907819, + 0.06861998972363427, 0.23802538347833696, 0.19769987611740356}; + +} // namespace + +TEST_CASE("BlackBoxPrimitiveEstimator") { + pluginplay::ModuleManager mm; + integrals::load_modules(mm); + auto& mod = mm.at("Black Box Primitive Pair Estimator"); + auto h2_basis = test::h2_sto3g_basis_set(); + + SECTION("H2 STO-3G") { + const auto& Kij = mod.run_as(h2_basis, h2_basis); + tensorwrapper::shape::Smooth shape{6, 6}; + tensorwrapper::buffer::Contiguous corr(h2_h2_corr, shape); + REQUIRE(Kij.buffer().approximately_equal(corr, 1e-8)); + } + + // Tests bra != ket + SECTION("Block of H2 STO-3G") { + simde::type::ao_basis_set h0_basis; + h0_basis.add_center(h2_basis.at(0)); + const auto& Kij = mod.run_as(h0_basis, h2_basis); + tensorwrapper::shape::Smooth shape{3, 6}; + std::vector buffer(h2_h2_corr.begin(), h2_h2_corr.begin() + 18); + tensorwrapper::buffer::Contiguous corr(buffer, shape); + REQUIRE(Kij.buffer().approximately_equal(corr, 1e-8)); + } +} \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp new file mode 100644 index 00000000..5a6c70e9 --- /dev/null +++ b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -0,0 +1,31 @@ +/* + * 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 "test_error.hpp" +#include + +using pt = integrals::property_types::PrimitivePairEstimator; +using namespace integrals::libint::test; +TEST_CASE("CauchySchwarzPrimitiveEstimator") { + pluginplay::ModuleManager mm; + integrals::load_modules(mm); + integrals::set_defaults(mm); + auto& mod = mm.at("CauchySchwarz Estimator"); + auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + + auto Q_ab = mod.run_as(bra0, bra1); + auto Q_cd = mod.run_as(ket0, ket1); + // std::cout << Q_ab << std::endl; +} \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp new file mode 100644 index 00000000..666503b1 --- /dev/null +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -0,0 +1,41 @@ +#include "test_error.hpp" +#include + +using eri4_pt = simde::ERI4; +using pt = integrals::property_types::Uncertainty; + +using namespace integrals::libint::test; + +TEST_CASE("PrimitiveErrorModel") { + pluginplay::ModuleManager mm; + integrals::load_modules(mm); + integrals::set_defaults(mm); + + auto& mod = mm.at("Primitive Error Model"); + + simde::type::v_ee_type v_ee{}; + using float_type = double; + SECTION("H2") { + auto h2_aos = test::h2_sto3g_basis_set(); + simde::type::aos_squared h2_aos2(h2_aos, h2_aos); + chemist::braket::BraKet mnls(h2_aos2, v_ee, h2_aos2); + auto error = mod.run_as(mnls, 1.0e-10); + + tensorwrapper::shape::Smooth shape({2, 2, 2, 2}); + std::vector data(shape.size()); + tensorwrapper::buffer::Contiguous buffer(std::move(data), shape); + REQUIRE(buffer.approximately_equal(error.buffer(), 1.0e-12)); + } + + SECTION("H2 2") { + float_type tol = 1.0e-10; + auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + + simde::type::aos_squared bra(bra0, bra1); + simde::type::aos_squared ket(ket0, ket1); + chemist::braket::BraKet mnls(bra, v_ee, ket); + auto error = mod.run_as(mnls, tol); + + auto corr = mm.at("Analytic Error").run_as(mnls, tol); + } +} \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/test_error.hpp b/tests/cxx/unit/integrals/libint/test_error.hpp new file mode 100644 index 00000000..f86d233b --- /dev/null +++ b/tests/cxx/unit/integrals/libint/test_error.hpp @@ -0,0 +1,77 @@ +#pragma once +#include "testing.hpp" + +namespace integrals::libint::test { + +/* After playing around, I found that if you set up a H2 dimer so that the + H2 molecules are 3.0 Angstroms apart, then the (03|12) integral has a + relatively large error (about 5.5E-9 using a tolerance of 1.0E-10). This + function prepares 4, 1-shell AO basis sets so that you can easily evaluate + that element. + */ +inline auto get_h2_dimer_0312_bases() { + simde::type::ao_basis_set bra0, bra1, ket0, ket1; + auto h2_2_aos = ::test::h2_sto3g_basis_set(); + + // Centers 0 and 1 stay put, center 2 is 0 translated, and 3 is 1 translated + bra0.add_center(h2_2_aos[0]); // 0 + bra1.add_center(h2_2_aos[1]); // 3 + ket0.add_center(h2_2_aos[1]); // 1 + ket1.add_center(h2_2_aos[0]); // 2 + + const double r = 3.0; + auto r0 = h2_2_aos[0].center(); + auto r1 = h2_2_aos[1].center(); + + bra1[0].center().x() = r1.x() + r; + bra1[0].center().y() = r1.y() + r; + bra1[0].center().z() = r1.z() + r; + ket1[0].center().x() = r0.x() + r; + ket1[0].center().y() = r0.y() + r; + ket1[0].center().z() = r0.z() + r; + return std::make_tuple(bra0, bra1, ket0, ket1); +} + +/* Given four 1 shell basis sets, this function will create the shell quartet + in the AO basis set by contracting the primitive integrals with the + primitive coefficients. This is a sanity check more than anything else. + */ +template +inline auto manual_contract_shell(const BasisType& bra0, const BasisType& bra1, + const BasisType& ket0, const BasisType& ket1, + pluginplay::ModuleManager& mm) { + auto dec_mod = mm.at("Decontract Basis Set"); + using dec_pt = integrals::property_types::DecontractBasisSet; + auto bra0_prims = dec_mod.run_as(bra0); + auto bra1_prims = dec_mod.run_as(bra1); + auto ket0_prims = dec_mod.run_as(ket0); + auto ket1_prims = dec_mod.run_as(ket1); + simde::type::v_ee_type v_ee{}; + simde::type::aos_squared bra_prims(bra0_prims, bra1_prims); + simde::type::aos_squared ket_prims(ket0_prims, ket1_prims); + chemist::braket::BraKet mnls(bra_prims, v_ee, ket_prims); + auto ints = mm.at("ERI4").run_as(mnls); + + const auto& buffer = tensorwrapper::buffer::make_contiguous(ints.buffer()); + + // TODO: expand to higher angular momentum + double value = 0.0; + for(std::size_t i = 0; i < bra0_prims.n_primitives(); ++i) { + auto ci = bra0.primitive(i).coefficient(); + for(std::size_t j = 0; j < bra1_prims.n_primitives(); ++j) { + auto cj = bra1.primitive(j).coefficient(); + for(std::size_t k = 0; k < ket0_prims.n_primitives(); ++k) { + auto ck = ket0.primitive(k).coefficient(); + for(std::size_t l = 0; l < ket1_prims.n_primitives(); ++l) { + auto cl = ket1.primitive(l).coefficient(); + auto erased = buffer.get_elem({i, j, k, l}); + auto val = wtf::fp::float_cast(erased); + value += ci * cj * ck * cl * val; + } + } + } + } + return value; +} + +} // namespace integrals::libint::test \ No newline at end of file diff --git a/tests/cxx/unit/integrals/testing.hpp b/tests/cxx/unit/integrals/testing.hpp index 812b05c3..b56a0c32 100644 --- a/tests/cxx/unit/integrals/testing.hpp +++ b/tests/cxx/unit/integrals/testing.hpp @@ -28,6 +28,12 @@ namespace test { +template +auto make_tw_buffer(FloatType value, tensorwrapper::shape::Smooth shape) { + std::vector data(1, value); + return tensorwrapper::buffer::Contiguous(std::move(data), shape); +} + template auto eigen_tensor_(const tensorwrapper::buffer::BufferBase& buffer, std::array extents, std::index_sequence) { @@ -183,6 +189,7 @@ inline simde::type::ao_basis_set h2_sto3g_basis_set() { using point_t = simde::type::point; using doubles_t = std::vector; + // Yes this is wrong, but it's what all the tests are expecting now auto mol = water_molecule(); point_t r0 = mol[0].as_nucleus(); point_t r1 = mol[1].as_nucleus(); diff --git a/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp b/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp new file mode 100644 index 00000000..d86f81c8 --- /dev/null +++ b/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp @@ -0,0 +1,99 @@ +/* + * 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 "../libint/test_error.hpp" +#include + +using namespace integrals::libint::test; +using namespace integrals::utils; + +using eri_pt = simde::ERI4; +using size_type = std::size_t; +using Catch::Approx; +using tensorwrapper::buffer::make_contiguous; + +/* Testing notes. + * + * All correct answers were generated by using a Python script. Since the + * script basically does the same thing as the C++ code, there is a chance we + * have replicated the same error in both places. However, since the Python + * script uses a different mechanism (i.e., Numpy) this is less likely. The + * script is: + * + * ``` + * import numpy as np + * eris = np.array(...) # Put ERI block here + * norm2 = np.einsum('mnls,mnls->', eris, eris) + * print(np.sqrt(norm2)) + * ``` + */ + +TEST_CASE("Rank2 Shell Norm") { + pluginplay::ModuleManager mm; + integrals::load_modules(mm); + integrals::set_defaults(mm); + + auto eri_mod = mm.at("ERI4"); + + auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + simde::type::aos_squared bra(bra0, bra1); + simde::type::aos_squared ket(ket0, ket1); + simde::type::v_ee_type v_ee{}; + + // Initially we test one quartet at a time + std::array offsets{0, 0, 0, 0}; + + SECTION("(ss|ss)") { + std::array naos{1, 1, 1, 1}; + chemist::braket::BraKet mnls(bra, v_ee, ket); + auto I = eri_mod.run_as(mnls); + auto& buffer = make_contiguous(I.buffer()); + auto norm = rank2_shell_norm(buffer, offsets, naos); + REQUIRE(norm == Approx(0.0002263495626484).epsilon(1.0e-8)); + } + + SECTION("(sp|ss)") { + bra.second.ao_basis_set().shell(0).l() = 1; + std::array naos{1, 3, 1, 1}; + chemist::braket::BraKet mnls(bra, v_ee, ket); + auto I = eri_mod.run_as(mnls); + auto& buffer = make_contiguous(I.buffer()); + auto norm = rank2_shell_norm(buffer, offsets, naos); + REQUIRE(norm == Approx(0.0007106429408912309).epsilon(1.0e-8)); + } + + SECTION("(sd|ps)") { + bra.second.ao_basis_set().shell(0).l() = 2; + ket.first.ao_basis_set().shell(0).l() = 1; + std::array naos{1, 6, 3, 1}; + chemist::braket::BraKet mnls(bra, v_ee, ket); + auto I = eri_mod.run_as(mnls); + auto& buffer = make_contiguous(I.buffer()); + auto norm = rank2_shell_norm(buffer, offsets, naos); + REQUIRE(norm == Approx(0.0021424444248565305).epsilon(1.0e-8)); + } + + SECTION("(sp|df)") { + bra.second.ao_basis_set().shell(0).l() = 1; + ket.first.ao_basis_set().shell(0).l() = 2; + ket.second.ao_basis_set().shell(0).l() = 3; + std::array naos{1, 3, 6, 10}; + chemist::braket::BraKet mnls(bra, v_ee, ket); + auto I = eri_mod.run_as(mnls); + auto& buffer = make_contiguous(I.buffer()); + auto norm = rank2_shell_norm(buffer, offsets, naos); + REQUIRE(norm == Approx(0.0035472410738408154).epsilon(1.0e-8)); + } +} \ No newline at end of file From b31e89a9e06bb4825bf1101c658592ca1d8ddfdd Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 10 Feb 2026 08:40:26 -0600 Subject: [PATCH 02/14] backukp --- src/integrals/libint/libint.cpp | 5 +- .../libint/primitive_error_model.cpp | 35 ++++++++--- .../libint/primitive_error_model.cpp | 13 ++++- .../cxx/unit/integrals/libint/test_error.hpp | 58 +++++++++++++++++-- 4 files changed, 96 insertions(+), 15 deletions(-) diff --git a/src/integrals/libint/libint.cpp b/src/integrals/libint/libint.cpp index 4395a1e5..b63618e8 100644 --- a/src/integrals/libint/libint.cpp +++ b/src/integrals/libint/libint.cpp @@ -192,8 +192,11 @@ EXTERN_LIBINT(aos_squared, v_ee_type, aos_squared); #undef EXTERN_LIBINT void set_defaults(pluginplay::ModuleManager& mm) { - mm.change_submod("Primitive Error Model", "Primitive Pair Estimator", + mm.change_submod("Primitive Error Model", + "Black Box Primitive Pair Estimator", "Black Box Primitive Pair Estimator"); + mm.change_submod("Primitive Error Model", "Primitive Pair Estimator", + "CauchySchwarz Estimator"); mm.change_submod("CauchySchwarz Estimator", "Decontract Basis Set", "Decontract Basis Set"); mm.change_submod("CauchySchwarz Estimator", "ERI4", "ERI4"); diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp index 5ac49c0a..3ffddf42 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -10,7 +10,8 @@ const auto desc = "Compute uncertainty estimates for ERI4 integrals"; // Sums contributions Q_abQ_cd when Q_ab or Q_cd is below tol template auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, -TensorType&& Q_ab, TensorType&& Q_cd, double tol){ +TensorType&& Q_ab, TensorType&& Q_cd, +TensorType&& K_ab, TensorType&& K_cd, double tol){ // Get the number of primitives in each shell const auto n_prim0 = shells[0].n_primitives(); const auto n_prim1 = shells[1].n_primitives(); @@ -37,9 +38,12 @@ TensorType&& Q_ab, TensorType&& Q_cd, double tol){ // Was the Q_ab element neglected for primitive pair (p0, p1)? std::vector i01{p0, p1}; + auto K_01 = float_cast(K_ab.get_elem(i01)); + auto K_01_mag = std::fabs(K_01); + bool K_01_neglected = K_01_mag < tol; auto Q_01 = float_cast(Q_ab.get_elem(i01)); auto Q_01_mag = std::fabs(Q_01); - + for(prim[2] = 0; prim[2] < n_prim2; ++prim[2]){ const auto p2 = prim_offsets[2] + prim[2]; @@ -48,12 +52,18 @@ TensorType&& Q_ab, TensorType&& Q_cd, double tol){ // Was the Q_cd element neglected for pair (p2, p3)? std::vector i23{p2, p3}; + auto K_23 = float_cast(K_cd.get_elem(i23)); + auto K_23_mag = std::fabs(K_23); + bool K_23_neglected = K_23_mag < tol; auto Q_23 = float_cast(Q_cd.get_elem(i23)); auto Q_23_mag = std::fabs(Q_23); - + + auto prod_neglected = K_01_mag * K_23_mag < tol; + // If either was neglected we pick up an error Q_01 * Q_23 - if(Q_01_mag * Q_23_mag < tol) + if(K_01_neglected || K_23_neglected || prod_neglected){ error += Q_01_mag * Q_23_mag; + } } // End loop over prim[3] } // End loop over prim[2] } // End loop over prim[1] @@ -108,7 +118,10 @@ MODULE_CTOR(PrimitiveErrorModel) { description(desc); // TODO citation for Chemist paper - add_submodule("Primitive Pair Estimator") + add_submodule("Black Box Primitive Pair Estimator") + .set_description("The module used to estimate the contributions of " + "primitive pairs to the overall integral values"); + add_submodule("Primitive Pair Estimator") .set_description("The module used to estimate the contributions of " "primitive pairs to the overall integral values"); } @@ -121,15 +134,23 @@ MODULE_RUN(PrimitiveErrorModel) { const auto ket0 = braket.ket().first.ao_basis_set(); const auto ket1 = braket.ket().second.ao_basis_set(); - // Get the Q_ab and Q_cd tensors + // Get the Q_ab and Q_cd tensors (used to estimate error) auto& estimator = submods.at("Primitive Pair Estimator"); const auto& Q_ab_tw = estimator.run_as(bra0, bra1); const auto& Q_cd_tw = estimator.run_as(ket0, ket1); + // Get the K_ab and K_cd tensors (used to determine which primitives are + // neglected) + auto& estimator2 = submods.at("Black Box Primitive Pair Estimator"); + const auto& K_ab_tw = estimator2.run_as(bra0, bra1); + const auto& K_cd_tw = estimator2.run_as(ket0, ket1); + // XXX: Workaround for needing the contiguous buffers to access elements using tensorwrapper::buffer::make_contiguous; const auto& Q_ab = make_contiguous(Q_ab_tw.buffer()); const auto& Q_cd = make_contiguous(Q_cd_tw.buffer()); + const auto& K_ab = make_contiguous(K_ab_tw.buffer()); + const auto& K_cd = make_contiguous(K_cd_tw.buffer()); // Get the number of AOs in each basis set const auto n_mu = bra0.n_aos(); @@ -174,7 +195,7 @@ MODULE_RUN(PrimitiveErrorModel) { shells[3] = ket1.shell(shell_i[3]); auto shell_error = - shell_block_error(shells, prim_offset, Q_ab, Q_cd, tol); + shell_block_error(shells, prim_offset, Q_ab, Q_cd, K_ab, K_cd, tol); fill_ao_block(shells, ao_offsets, buffer, shell_error); prim_offset[3] += shells[3].n_primitives(); diff --git a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp index 666503b1..42a45938 100644 --- a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -30,12 +30,23 @@ TEST_CASE("PrimitiveErrorModel") { SECTION("H2 2") { float_type tol = 1.0e-10; auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + auto value = manual_contract_shell(bra0, bra1, ket0, ket1, mm); + std::cout << "Screened value: " << value << std::endl; simde::type::aos_squared bra(bra0, bra1); simde::type::aos_squared ket(ket0, ket1); chemist::braket::BraKet mnls(bra, v_ee, ket); + auto error = mod.run_as(mnls, tol); + auto corr = mm.at("Analytic Error").run_as(mnls, tol); + + auto eri4_mod = mm.at("ERI4").unlocked_copy(); + eri4_mod.change_input("Threshold", float_type(1E-16)); + auto i_value = eri4_mod.run_as(mnls); + std::cout << "Corr Screened Value: " << i_value << std::endl; - auto corr = mm.at("Analytic Error").run_as(mnls, tol); + // Our value: 0.0002263495591894 + // Libint's value: 0.0002263440759384 + // No screening: 0.0002263495626484 } } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/test_error.hpp b/tests/cxx/unit/integrals/libint/test_error.hpp index f86d233b..32cfa4e1 100644 --- a/tests/cxx/unit/integrals/libint/test_error.hpp +++ b/tests/cxx/unit/integrals/libint/test_error.hpp @@ -32,6 +32,20 @@ inline auto get_h2_dimer_0312_bases() { return std::make_tuple(bra0, bra1, ket0, ket1); } +template +double black_box_pair_metric(const T& prim_i, const T& prim_j) { + auto ci = prim_i.coefficient(); + auto cj = prim_j.coefficient(); + auto zeta_i = prim_i.exponent(); + auto zeta_j = prim_j.exponent(); + auto dx = prim_i.center().x() - prim_j.center().x(); + auto dy = prim_i.center().y() - prim_j.center().y(); + auto dz = prim_i.center().z() - prim_j.center().z(); + auto dr = std::sqrt(dx * dx + dy * dy + dz * dz); + auto exponent = (-zeta_i * zeta_j / (zeta_i + zeta_j)) * dr * dr; + return ci * cj * std::exp(exponent); +} + /* Given four 1 shell basis sets, this function will create the shell quartet in the AO basis set by contracting the primitive integrals with the primitive coefficients. This is a sanity check more than anything else. @@ -51,27 +65,59 @@ inline auto manual_contract_shell(const BasisType& bra0, const BasisType& bra1, simde::type::aos_squared ket_prims(ket0_prims, ket1_prims); chemist::braket::BraKet mnls(bra_prims, v_ee, ket_prims); auto ints = mm.at("ERI4").run_as(mnls); + // auto K_mod = mm.at("Black Box Primitive Pair Estimator"); + // using prim_pt = integrals::property_types::PrimitivePairEstimator; + // auto K_ab = K_mod.run_as(bra0, bra1); + // auto K_cd = K_mod.run_as(ket0, ket1); const auto& buffer = tensorwrapper::buffer::make_contiguous(ints.buffer()); + // const auto& K01 = + // tensorwrapper::buffer::make_contiguous(K_ab.buffer()); const auto& K23 = + // tensorwrapper::buffer::make_contiguous(K_cd.buffer()); // TODO: expand to higher angular momentum double value = 0.0; + double tol = 1e-10; + double error = 0.0; + using wtf::fp::float_cast; for(std::size_t i = 0; i < bra0_prims.n_primitives(); ++i) { - auto ci = bra0.primitive(i).coefficient(); + const auto prim_i = bra0.primitive(i); + auto ci = prim_i.coefficient(); + for(std::size_t j = 0; j < bra1_prims.n_primitives(); ++j) { - auto cj = bra1.primitive(j).coefficient(); + const auto prim_j = bra1.primitive(j); + auto cj = prim_j.coefficient(); + + auto kij = black_box_pair_metric(prim_i, prim_j); + + bool is01good = kij > tol; + for(std::size_t k = 0; k < ket0_prims.n_primitives(); ++k) { - auto ck = ket0.primitive(k).coefficient(); + const auto prim_k = ket0.primitive(k); + auto ck = prim_k.coefficient(); + for(std::size_t l = 0; l < ket1_prims.n_primitives(); ++l) { - auto cl = ket1.primitive(l).coefficient(); + const auto prim_l = ket1.primitive(l); + auto cl = prim_l.coefficient(); + + auto kkl = black_box_pair_metric(prim_k, prim_l); + bool is23good = kkl > tol; + bool both_good = (kij * kkl) > tol; + auto erased = buffer.get_elem({i, j, k, l}); auto val = wtf::fp::float_cast(erased); - value += ci * cj * ck * cl * val; + double Iijk = ci * cj * ck * cl * val; + if(is01good && is23good && both_good) { + value += Iijk; + } else { + error += std::fabs(Iijk); + } } } } } - return value; + std::cout << "Computed error: " << simde::type::tensor(error) << std::endl; + return simde::type::tensor(value); } } // namespace integrals::libint::test \ No newline at end of file From e7841d6089dec393e5f53e335e2c71fd96a7c6ad Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Thu, 19 Feb 2026 21:51:05 -0600 Subject: [PATCH 03/14] backup --- .../libint/black_box_primitive_estimator.cpp | 93 ++++++++++++------- 1 file changed, 62 insertions(+), 31 deletions(-) diff --git a/src/integrals/libint/black_box_primitive_estimator.cpp b/src/integrals/libint/black_box_primitive_estimator.cpp index 4a72119c..15b56def 100644 --- a/src/integrals/libint/black_box_primitive_estimator.cpp +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -36,43 +36,74 @@ MODULE_CTOR(BlackBoxPrimitiveEstimator) { MODULE_RUN(BlackBoxPrimitiveEstimator) { const auto&& [bra_basis, ket_basis] = pt::unwrap_inputs(inputs); - const auto n_bra = bra_basis.n_primitives(); - const auto n_ket = ket_basis.n_primitives(); + const auto n_shells_bra = bra_basis.n_shells(); + const auto n_prims_bra = bra_basis.n_primitives(); + const auto n_shells_ket = ket_basis.n_shells(); + const auto n_prims_ket = ket_basis.n_primitives(); + + // TODO: Our basis really needs to handle normalization better... + auto bra = make_libint_basis_set(bra_basis); + auto ket = make_libint_basis_set(ket_basis); // TODO: Use floating point type of the basis sets using float_type = double; - std::vector buffer(n_bra * n_ket, 0.0); + std::vector buffer(n_prims_bra * n_prims_ket, 0.0); using iter_type = std::decay_t; // Type of the loop index - - for(iter_type bra_i = 0; bra_i < n_bra; ++bra_i) { - const auto bra_prim_i = bra_basis.primitive(bra_i); - const auto bra_zeta = bra_prim_i.exponent(); - const auto bra_coeff = std::fabs(bra_prim_i.coefficient()); - const auto bra_center = bra_prim_i.center(); - const auto bra_offset = bra_i * n_ket; - - for(iter_type ket_i = 0; ket_i < n_ket; ++ket_i) { - const auto ket_prim_i = ket_basis.primitive(ket_i); - const auto ket_zeta = ket_prim_i.exponent(); - const auto ket_coeff = std::fabs(ket_prim_i.coefficient()); - const auto ket_center = ket_prim_i.center(); - - // This is "K bar" in Eq. 11 in the SI of the Chemist paper - const auto dist = (bra_center - ket_center).magnitude(); - const auto dist2 = dist * dist; - const auto ratio = bra_zeta * ket_zeta / (bra_zeta + ket_zeta); - const auto coeff = bra_coeff * ket_coeff; - buffer[bra_offset + ket_i] = coeff * std::exp(-1.0 * ratio * dist2); + iter_type bra_counter = 0; + + for(iter_type bra_shell_i = 0; bra_shell_i < n_shells_bra; ++bra_shell_i) { + const auto& bra_shell = bra.at(bra_shell_i); + assert(bra_shell.contr.size() == 1); + const auto& bra_coeff = bra_shell.contr[0].coeff; + const auto& bra_alpha = bra_shell.contr[0].alpha; + const auto n_prims_bra_shell = bra_coeff.size(); + + for(iter_type bra_prim_i = 0; bra_prim_i < n_prims_bra_shell; + ++bra_prim_i) { + const auto bra_zeta = bra_alpha[bra_prim_i]; + const auto bra_coeff = std::fabs(bra_coeff[bra_prim_i]); + const auto bra_center = bra_shell.O; + const auto bra_offset = (bra_counter + bra_prim_i) * n_prims_ket; + + iter_type ket_counter = 0; + for(iter_type ket_shell_i = 0; ket_shell_i < n_shells_ket; + ++ket_shell_i) { + const auto& ket_shell = ket.at(ket_shell_i); + assert(ket_shell.contr.size() == 1); + const auto& ket_coeff = ket_shell.contr[0].coeff; + const auto& ket_alpha = ket_shell.contr[0].alpha; + const auto n_prims_ket_shell = ket_coeff.size(); + + for(iter_type ket_prim_i = 0; ket_prim_i < n_prims_ket_shell; + ++ket_prim_i) { + const auto ket_zeta = ket_alpha[ket_prim_i]; + const auto ket_coeff = std::fabs(ket_coeff[ket_prim_i]); + const auto ket_center = ket_shell.O; + const auto ket_offset = ket_counter + ket_prim_i; + + // This is "K bar" in Eq. 11 in the SI of the Chemist paper + const auto dx = bra_center.x - ket_center.x; + const auto dy = bra_center.y - ket_center.y; + const auto dz = bra_center.z - ket_center.z; + const auto dr2 = dx * dx + dy * dy + dz * dz; + + const auto ratio = + bra_zeta * ket_zeta / (bra_zeta + ket_zeta); + const auto coeff = bra_coeff * ket_coeff; + buffer[bra_offset + ket_offset] = + coeff * std::exp(-1.0 * ratio * dist2); + } + ket_counter += n_prims_ket_shell; + } + bra_counter += n_prims_bra_shell; } - } + tensorwrapper::shape::Smooth shape({n_bra, n_ket}); + tensorwrapper::buffer::Contiguous tw_buffer(std::move(buffer), shape); + simde::type::tensor rv(shape, std::move(tw_buffer)); - tensorwrapper::shape::Smooth shape({n_bra, n_ket}); - tensorwrapper::buffer::Contiguous tw_buffer(std::move(buffer), shape); - simde::type::tensor rv(shape, std::move(tw_buffer)); - - auto result = results(); - return pt::wrap_results(result, rv); -} + auto result = results(); + return pt::wrap_results(result, rv); + } } // namespace integrals::libint \ No newline at end of file From 2b2732a9779b078fab276e0d8714cbe615c253bd Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 20 Feb 2026 15:56:52 -0600 Subject: [PATCH 04/14] refactor estimator --- .../libint/black_box_primitive_estimator.cpp | 138 +++++++++++------- .../libint/black_box_primitive_estimator.cpp | 28 ++-- 2 files changed, 101 insertions(+), 65 deletions(-) diff --git a/src/integrals/libint/black_box_primitive_estimator.cpp b/src/integrals/libint/black_box_primitive_estimator.cpp index 15b56def..38b1c50f 100644 --- a/src/integrals/libint/black_box_primitive_estimator.cpp +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -21,10 +21,41 @@ namespace integrals::libint { namespace { -const auto desc = ""; +const auto desc = R"( +Libint Black Box Primitive Pair Estimator +========================================= +This module computes the matrix : math :`K_{ij}` where: + +.. math:: + + K_{ij} = c_i c_j \exp\left(-\frac{\zeta_i \zeta_j}{\zeta_i + \zeta_j} + |\mathbf{R}_i - \mathbf{R}_j|^2\right) + +This is how Libint2 estimates the contribution of a pair of primitives to +an integral. + +)"; + +// Computes square of the distance between points a and b +// (T should be libint::Atom like) +template +auto distance_squared(T&& a, T&& 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; } +template +auto compute_k(T zeta_i, T zeta_j, T coeff_i, T coeff_j, T dr2) { + const auto num = -zeta_i * zeta_j; + const auto denom = zeta_i + zeta_j; + const auto ratio = num / denom; + return coeff_i * coeff_j * std::exp(ratio * dr2); +} + +} // namespace using pt = integrals::property_types::PrimitivePairEstimator; MODULE_CTOR(BlackBoxPrimitiveEstimator) { @@ -34,76 +65,81 @@ MODULE_CTOR(BlackBoxPrimitiveEstimator) { } MODULE_RUN(BlackBoxPrimitiveEstimator) { - const auto&& [bra_basis, ket_basis] = pt::unwrap_inputs(inputs); + const auto&& [bra, ket] = pt::unwrap_inputs(inputs); - const auto n_shells_bra = bra_basis.n_shells(); - const auto n_prims_bra = bra_basis.n_primitives(); - const auto n_shells_ket = ket_basis.n_shells(); - const auto n_prims_ket = ket_basis.n_primitives(); + using iter_type = std::size_t; // Type of each loop index + using index_array = std::array; // Type of a pair of indices + using index_vector = std::vector; // Type of a vector of indices + using float_type = double; // TODO: Get from basis sets - // TODO: Our basis really needs to handle normalization better... - auto bra = make_libint_basis_set(bra_basis); - auto ket = make_libint_basis_set(ket_basis); + index_array shells{0, 0}; // shells[0]/shells[1] indexes bra/ket shell + index_array n_shells{bra.n_shells(), ket_.n_shells()} // Number of shells - // TODO: Use floating point type of the basis sets - using float_type = double; - std::vector buffer(n_prims_bra * n_prims_ket, 0.0); + index_array prims{0, 0}; // prims[0]/prims[1] indexes bra/ket primitive + index_array n_prims{bra.n_primitives(), ket.n_primitives()}; + index_array offsets{0, 0}; // Offset to the first primitive of the shell + index_vector abs_prim{0, 0}; // Absolute indices of the primitives - using iter_type = std::decay_t; // Type of the loop index - iter_type bra_counter = 0; + // Will be the result + tensorwrapper::shape::Smooth shape({n_prims[0], n_prims[1]}); + std::vector data(shape.size(), 0.0); + tensorwrapper::buffer::Contiguous buffer(std::move(data), shape); - for(iter_type bra_shell_i = 0; bra_shell_i < n_shells_bra; ++bra_shell_i) { - const auto& bra_shell = bra.at(bra_shell_i); + // For now use the libint basis sets because they're properly normalized + // TODO: Our basis really needs to handle normalization better... + auto bra_libint = make_libint_basis_set(bra); + auto ket_libint = make_libint_basis_set(ket); + + for(shells[0] = 0; shells[0] < n_shells[0]; ++shells[0]) { + const auto& bra_shell = bra_libint.at(shells[0]); assert(bra_shell.contr.size() == 1); const auto& bra_coeff = bra_shell.contr[0].coeff; const auto& bra_alpha = bra_shell.contr[0].alpha; const auto n_prims_bra_shell = bra_coeff.size(); - for(iter_type bra_prim_i = 0; bra_prim_i < n_prims_bra_shell; - ++bra_prim_i) { - const auto bra_zeta = bra_alpha[bra_prim_i]; - const auto bra_coeff = std::fabs(bra_coeff[bra_prim_i]); + for(prims[0] = 0; prims[0] < n_prims_bra_shell; ++prims[0]) { + const auto zeta0 = bra_alpha[prims[0]]; + const auto coeff0 = std::fabs(bra_coeff[prims[0]]); const auto bra_center = bra_shell.O; - const auto bra_offset = (bra_counter + bra_prim_i) * n_prims_ket; + abs_prim[0] = offsets[0] + prims[0]; - iter_type ket_counter = 0; - for(iter_type ket_shell_i = 0; ket_shell_i < n_shells_ket; - ++ket_shell_i) { - const auto& ket_shell = ket.at(ket_shell_i); + offsets[1] = 0; + for(shells[1] = 0; shells[1] < n_shells[1]; ++shells[1]) { + const auto& ket_shell = ket.at(shells[1]); assert(ket_shell.contr.size() == 1); const auto& ket_coeff = ket_shell.contr[0].coeff; const auto& ket_alpha = ket_shell.contr[0].alpha; const auto n_prims_ket_shell = ket_coeff.size(); - for(iter_type ket_prim_i = 0; ket_prim_i < n_prims_ket_shell; - ++ket_prim_i) { - const auto ket_zeta = ket_alpha[ket_prim_i]; - const auto ket_coeff = std::fabs(ket_coeff[ket_prim_i]); + for(prims[1] = 0; prims[1] < n_prims_ket_shell; ++prims[1]) { + const auto zeta1 = ket_alpha[prims[1]]; + const auto coeff1 = std::fabs(ket_coeff[prims[1]]); const auto ket_center = ket_shell.O; - const auto ket_offset = ket_counter + ket_prim_i; + abs_prim[1] = ket_counter + ket_prim_i; // This is "K bar" in Eq. 11 in the SI of the Chemist paper - const auto dx = bra_center.x - ket_center.x; - const auto dy = bra_center.y - ket_center.y; - const auto dz = bra_center.z - ket_center.z; - const auto dr2 = dx * dx + dy * dy + dz * dz; - - const auto ratio = - bra_zeta * ket_zeta / (bra_zeta + ket_zeta); - const auto coeff = bra_coeff * ket_coeff; - buffer[bra_offset + ket_offset] = - coeff * std::exp(-1.0 * ratio * dist2); - } + auto dr2 = distance_squared(bra_center, ket_center); + auto k01 = compute_k(zeta0, zeta1, coeff0, coeff1, dr2); + buffer.set_element(abs_prim, k01); + } // loop over ket primitives + ket_counter += n_prims_ket_shell; - } - bra_counter += n_prims_bra_shell; - } - tensorwrapper::shape::Smooth shape({n_bra, n_ket}); - tensorwrapper::buffer::Contiguous tw_buffer(std::move(buffer), shape); - simde::type::tensor rv(shape, std::move(tw_buffer)); - - auto result = results(); - return pt::wrap_results(result, rv); - } + } // loop over ket shells + + // We ran over all ket primitives, so counter should be done too + assert(ket_counter == n_prims[1]); + } // loop over bra primitives + + bra_counter += n_prims_bra_shell; + } // loop over bra shells + + // We ran over all bra primitives, so counter should be done too + assert(bra_counter == n_prims[0]); + + simde::type::tensor rv(shape, std::move(buffer)); + + auto result = results(); + return pt::wrap_results(result, rv); +} } // namespace integrals::libint \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp index de2a07f6..65df1b71 100644 --- a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp +++ b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp @@ -20,21 +20,21 @@ using pt = integrals::property_types::PrimitivePairEstimator; namespace { -std::vector h2_h2_corr{ - 0.023817430147884473, 0.08261663936778646, 0.06861998972363427, - 1.4555047643750448e-05, 0.00844595180076824, 0.03423471319097009, - 0.08261663936778646, 0.28657621993836907, 0.23802538347833696, - 0.00844595180076824, 0.07444365105885908, 0.13404298325096972, - 0.06861998972363427, 0.23802538347833696, 0.19769987611740356, - 0.03423471319097009, 0.13404298325096972, 0.1372685140907819, - 1.4555047643750448e-05, 0.00844595180076824, 0.03423471319097009, - 0.023817430147884473, 0.08261663936778646, 0.06861998972363427, - 0.00844595180076824, 0.07444365105885908, 0.13404298325096972, - 0.08261663936778646, 0.28657621993836907, 0.23802538347833696, - 0.03423471319097009, 0.13404298325096972, 0.1372685140907819, - 0.06861998972363427, 0.23802538347833696, 0.19769987611740356}; +// std::vector h2_h2_corr{ +// 0.023817430147884473, 0.08261663936778646, 0.06861998972363427, +// 1.4555047643750448e-05, 0.00844595180076824, 0.03423471319097009, +// 0.08261663936778646, 0.28657621993836907, 0.23802538347833696, +// 0.00844595180076824, 0.07444365105885908, 0.13404298325096972, +// 0.06861998972363427, 0.23802538347833696, 0.19769987611740356, +// 0.03423471319097009, 0.13404298325096972, 0.1372685140907819, +// 1.4555047643750448e-05, 0.00844595180076824, 0.03423471319097009, +// 0.023817430147884473, 0.08261663936778646, 0.06861998972363427, +// 0.00844595180076824, 0.07444365105885908, 0.13404298325096972, +// 0.08261663936778646, 0.28657621993836907, 0.23802538347833696, +// 0.03423471319097009, 0.13404298325096972, 0.1372685140907819, +// 0.06861998972363427, 0.23802538347833696, 0.19769987611740356}; -} // namespace +// } // namespace TEST_CASE("BlackBoxPrimitiveEstimator") { pluginplay::ModuleManager mm; From ff1ee5bdc025dbd93fe13f86b112adbb2ad2e77c Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Sun, 22 Feb 2026 11:40:22 -0600 Subject: [PATCH 05/14] uses normalized coeffs, refactor testing infrastructure --- .../libint/black_box_primitive_estimator.cpp | 38 +++-- .../libint/primitive_error_model.cpp | 111 +++++++------- .../ao_integrals/ao_integrals_driver.cpp | 13 +- .../integrals/ao_integrals/coulomb_metric.cpp | 14 +- .../integrals/ao_integrals/df_integral.cpp | 14 +- .../ao_integrals/j_density_fitted.cpp | 16 +- .../integrals/ao_integrals/j_four_center.cpp | 16 +- .../ao_integrals/k_density_fitted.cpp | 16 +- .../integrals/ao_integrals/k_four_center.cpp | 16 +- .../integrals/ao_integrals/test_uq_driver.cpp | 8 +- .../unit/integrals/libint/analytic_error.cpp | 30 ++-- .../libint/black_box_primitive_estimator.cpp | 75 ++++----- .../cauchy_schwarz_primitive_estimator.cpp | 56 ++++--- .../detail_/test_make_libint_basis_set.cpp | 6 +- .../libint/primitive_error_model.cpp | 102 ++++++------ .../libint/test_arbitrary_operator.cpp | 47 +++--- tests/cxx/unit/integrals/libint/test_eri2.cpp | 17 +- tests/cxx/unit/integrals/libint/test_eri3.cpp | 17 +- tests/cxx/unit/integrals/libint/test_eri4.cpp | 17 +- .../cxx/unit/integrals/libint/test_error.hpp | 123 --------------- .../unit/integrals/libint/test_kinetic.cpp | 17 +- .../unit/integrals/libint/test_nuclear.cpp | 16 +- .../unit/integrals/libint/test_overlap.cpp | 17 +- .../{testing.hpp => testing/ao_bases.hpp} | 145 +++--------------- .../cxx/unit/integrals/testing/molecules.hpp | 25 +++ .../unit/integrals/testing/shell_quartets.hpp | 36 +++++ tests/cxx/unit/integrals/testing/testing.hpp | 82 ++++++++++ .../unit/integrals/utils/rank2_shell_norm.cpp | 4 +- .../utils/test_decontract_basis_set.cpp | 11 +- 29 files changed, 518 insertions(+), 587 deletions(-) delete mode 100644 tests/cxx/unit/integrals/libint/test_error.hpp rename tests/cxx/unit/integrals/{testing.hpp => testing/ao_bases.hpp} (54%) create mode 100644 tests/cxx/unit/integrals/testing/molecules.hpp create mode 100644 tests/cxx/unit/integrals/testing/shell_quartets.hpp create mode 100644 tests/cxx/unit/integrals/testing/testing.hpp diff --git a/src/integrals/libint/black_box_primitive_estimator.cpp b/src/integrals/libint/black_box_primitive_estimator.cpp index 38b1c50f..e2b85567 100644 --- a/src/integrals/libint/black_box_primitive_estimator.cpp +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -14,6 +14,7 @@ * limitations under the License. */ +#include "detail_/make_libint_basis_set.hpp" #include "libint.hpp" #include #include @@ -35,15 +36,18 @@ This module computes the matrix : math :`K_{ij}` where: This is how Libint2 estimates the contribution of a pair of primitives to an integral. +N.B. The algorithm assumes that the bra and ket are different. If they are the +same, we can save time by using the fact that the matrix is symmetric. + )"; // Computes square of the distance between points a and b // (T should be libint::Atom like) template auto distance_squared(T&& a, T&& b) { - auto dx = a.x - b.x; - auto dy = a.y - b.y; - auto dz = a.z - b.z; + auto dx = a[0] - b[0]; + auto dy = a[1] - b[1]; + auto dz = a[2] - b[2]; return dx * dx + dy * dy + dz * dz; } @@ -73,7 +77,7 @@ MODULE_RUN(BlackBoxPrimitiveEstimator) { using float_type = double; // TODO: Get from basis sets index_array shells{0, 0}; // shells[0]/shells[1] indexes bra/ket shell - index_array n_shells{bra.n_shells(), ket_.n_shells()} // Number of shells + index_array n_shells{bra.n_shells(), ket.n_shells()}; // Number of shells index_array prims{0, 0}; // prims[0]/prims[1] indexes bra/ket primitive index_array n_prims{bra.n_primitives(), ket.n_primitives()}; @@ -87,14 +91,14 @@ MODULE_RUN(BlackBoxPrimitiveEstimator) { // For now use the libint basis sets because they're properly normalized // TODO: Our basis really needs to handle normalization better... - auto bra_libint = make_libint_basis_set(bra); - auto ket_libint = make_libint_basis_set(ket); + auto bra_libint = detail_::make_libint_basis_set(bra); + auto ket_libint = detail_::make_libint_basis_set(ket); for(shells[0] = 0; shells[0] < n_shells[0]; ++shells[0]) { const auto& bra_shell = bra_libint.at(shells[0]); - assert(bra_shell.contr.size() == 1); + assert(bra_shell.contr.size() == 1); // No general contraction support const auto& bra_coeff = bra_shell.contr[0].coeff; - const auto& bra_alpha = bra_shell.contr[0].alpha; + const auto& bra_alpha = bra_shell.alpha; const auto n_prims_bra_shell = bra_coeff.size(); for(prims[0] = 0; prims[0] < n_prims_bra_shell; ++prims[0]) { @@ -105,36 +109,36 @@ MODULE_RUN(BlackBoxPrimitiveEstimator) { offsets[1] = 0; for(shells[1] = 0; shells[1] < n_shells[1]; ++shells[1]) { - const auto& ket_shell = ket.at(shells[1]); - assert(ket_shell.contr.size() == 1); + const auto& ket_shell = ket_libint.at(shells[1]); + assert(ket_shell.contr.size() == 1); // No general contractions const auto& ket_coeff = ket_shell.contr[0].coeff; - const auto& ket_alpha = ket_shell.contr[0].alpha; + const auto& ket_alpha = ket_shell.alpha; const auto n_prims_ket_shell = ket_coeff.size(); for(prims[1] = 0; prims[1] < n_prims_ket_shell; ++prims[1]) { const auto zeta1 = ket_alpha[prims[1]]; const auto coeff1 = std::fabs(ket_coeff[prims[1]]); const auto ket_center = ket_shell.O; - abs_prim[1] = ket_counter + ket_prim_i; + abs_prim[1] = offsets[1] + prims[1]; // This is "K bar" in Eq. 11 in the SI of the Chemist paper auto dr2 = distance_squared(bra_center, ket_center); auto k01 = compute_k(zeta0, zeta1, coeff0, coeff1, dr2); - buffer.set_element(abs_prim, k01); + buffer.set_elem(abs_prim, k01); } // loop over ket primitives - ket_counter += n_prims_ket_shell; + offsets[1] += n_prims_ket_shell; } // loop over ket shells // We ran over all ket primitives, so counter should be done too - assert(ket_counter == n_prims[1]); + assert(offsets[1] == n_prims[1]); } // loop over bra primitives - bra_counter += n_prims_bra_shell; + offsets[0] += n_prims_bra_shell; } // loop over bra shells // We ran over all bra primitives, so counter should be done too - assert(bra_counter == n_prims[0]); + assert(offsets[0] == n_prims[0]); simde::type::tensor rv(shape, std::move(buffer)); diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp index 3ffddf42..873c11e1 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -1,6 +1,6 @@ #include "libint.hpp" -#include #include +#include namespace integrals::libint { namespace { @@ -9,9 +9,9 @@ const auto desc = "Compute uncertainty estimates for ERI4 integrals"; // Sums contributions Q_abQ_cd when Q_ab or Q_cd is below tol template -auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, -TensorType&& Q_ab, TensorType&& Q_cd, -TensorType&& K_ab, TensorType&& K_cd, double tol){ +auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, + TensorType&& Q_ab, TensorType&& Q_cd, TensorType&& K_ab, + TensorType&& K_cd, double tol) { // Get the number of primitives in each shell const auto n_prim0 = shells[0].n_primitives(); const auto n_prim1 = shells[1].n_primitives(); @@ -20,48 +20,47 @@ TensorType&& K_ab, TensorType&& K_cd, double tol){ using float_type = double; // TODO: get from Q_ab or Q_cd using iter_type = std::decay_t; // Type of the loop index - + // Create an array to hold the primitive indices std::array prim{0, 0, 0, 0}; - + // This is the accumulated error for this shell block double error = 0.0; - using wtf::fp::float_cast; - for(prim[0] = 0; prim[0] < n_prim0; ++prim[0]){ + for(prim[0] = 0; prim[0] < n_prim0; ++prim[0]) { const auto p0 = prim_offsets[0] + prim[0]; - for(prim[1] = 0; prim[1] < n_prim1; ++prim[1]){ + for(prim[1] = 0; prim[1] < n_prim1; ++prim[1]) { const auto p1 = prim_offsets[1] + prim[1]; - + // Was the Q_ab element neglected for primitive pair (p0, p1)? std::vector i01{p0, p1}; - auto K_01 = float_cast(K_ab.get_elem(i01)); - auto K_01_mag = std::fabs(K_01); + auto K_01 = float_cast(K_ab.get_elem(i01)); + auto K_01_mag = std::fabs(K_01); bool K_01_neglected = K_01_mag < tol; - auto Q_01 = float_cast(Q_ab.get_elem(i01)); - auto Q_01_mag = std::fabs(Q_01); + auto Q_01 = float_cast(Q_ab.get_elem(i01)); + auto Q_01_mag = std::fabs(Q_01); - for(prim[2] = 0; prim[2] < n_prim2; ++prim[2]){ + for(prim[2] = 0; prim[2] < n_prim2; ++prim[2]) { const auto p2 = prim_offsets[2] + prim[2]; - - for(prim[3] = 0; prim[3] < n_prim3; ++prim[3]){ + + for(prim[3] = 0; prim[3] < n_prim3; ++prim[3]) { const auto p3 = prim_offsets[3] + prim[3]; - + // Was the Q_cd element neglected for pair (p2, p3)? std::vector i23{p2, p3}; - auto K_23 = float_cast(K_cd.get_elem(i23)); + auto K_23 = float_cast(K_cd.get_elem(i23)); auto K_23_mag = std::fabs(K_23); bool K_23_neglected = K_23_mag < tol; - auto Q_23 = float_cast(Q_cd.get_elem(i23)); + auto Q_23 = float_cast(Q_cd.get_elem(i23)); auto Q_23_mag = std::fabs(Q_23); auto prod_neglected = K_01_mag * K_23_mag < tol; // If either was neglected we pick up an error Q_01 * Q_23 - if(K_01_neglected || K_23_neglected || prod_neglected){ + if(K_01_neglected || K_23_neglected || prod_neglected) { error += Q_01_mag * Q_23_mag; } } // End loop over prim[3] @@ -74,9 +73,8 @@ TensorType&& K_ab, TensorType&& K_cd, double tol){ // Sets all AO integrals in the shell block to the computed error template -void fill_ao_block(ShellsType&& shells, OffsetTypes&& ao_offsets, TensorType&& t, - double error){ - +void fill_ao_block(ShellsType&& shells, OffsetTypes&& ao_offsets, + TensorType&& t, double error) { // Get the number of AOs in each shell const auto n_aos0 = shells[0].size(); const auto n_aos1 = shells[1].size(); @@ -86,17 +84,17 @@ void fill_ao_block(ShellsType&& shells, OffsetTypes&& ao_offsets, TensorType&& t // Make an array to hold the AO indices using iter_type = std::decay_t; // Type of the loop index std::array ao{0, 0, 0, 0}; - - for(ao[0] = 0; ao[0] < n_aos0; ++ao[0]){ + + for(ao[0] = 0; ao[0] < n_aos0; ++ao[0]) { const auto a0 = ao_offsets[0] + ao[0]; - for(ao[1] = 0; ao[1] < n_aos1; ++ao[1]){ + for(ao[1] = 0; ao[1] < n_aos1; ++ao[1]) { const auto a1 = ao_offsets[1] + ao[1]; - - for(ao[2] = 0; ao[2] < n_aos2; ++ao[2]){ + + for(ao[2] = 0; ao[2] < n_aos2; ++ao[2]) { const auto a2 = ao_offsets[2] + ao[2]; - - for(ao[3] = 0; ao[3] < n_aos3; ++ao[3]){ + + for(ao[3] = 0; ao[3] < n_aos3; ++ao[3]) { const auto a3 = ao_offsets[3] + ao[3]; std::vector index{a0, a1, a2, a3}; @@ -121,7 +119,7 @@ MODULE_CTOR(PrimitiveErrorModel) { add_submodule("Black Box Primitive Pair Estimator") .set_description("The module used to estimate the contributions of " "primitive pairs to the overall integral values"); - add_submodule("Primitive Pair Estimator") + add_submodule("Primitive Pair Estimator") .set_description("The module used to estimate the contributions of " "primitive pairs to the overall integral values"); } @@ -135,13 +133,13 @@ MODULE_RUN(PrimitiveErrorModel) { const auto ket1 = braket.ket().second.ao_basis_set(); // Get the Q_ab and Q_cd tensors (used to estimate error) - auto& estimator = submods.at("Primitive Pair Estimator"); + auto& estimator = submods.at("Primitive Pair Estimator"); const auto& Q_ab_tw = estimator.run_as(bra0, bra1); const auto& Q_cd_tw = estimator.run_as(ket0, ket1); // Get the K_ab and K_cd tensors (used to determine which primitives are // neglected) - auto& estimator2 = submods.at("Black Box Primitive Pair Estimator"); + auto& estimator2 = submods.at("Black Box Primitive Pair Estimator"); const auto& K_ab_tw = estimator2.run_as(bra0, bra1); const auto& K_cd_tw = estimator2.run_as(ket0, ket1); @@ -153,10 +151,10 @@ MODULE_RUN(PrimitiveErrorModel) { const auto& K_cd = make_contiguous(K_cd_tw.buffer()); // Get the number of AOs in each basis set - const auto n_mu = bra0.n_aos(); - const auto n_nu = bra1.n_aos(); + const auto n_mu = bra0.n_aos(); + const auto n_nu = bra1.n_aos(); const auto n_lambda = ket0.n_aos(); - const auto n_sigma = ket1.n_aos(); + const auto n_sigma = ket1.n_aos(); // Make the return buffer using float_type = double; // TODO: get from Q_ab or Q_cd @@ -165,10 +163,10 @@ MODULE_RUN(PrimitiveErrorModel) { tensorwrapper::buffer::Contiguous buffer(std::move(raw_buffer), shape); // Make arrays to hold the indices and offsets - using iter_type = std::decay_t; // Type of the loop index + using iter_type = std::decay_t; // Type of the loop index using shell_view = decltype(bra0.shell(0)); - std::array n_shells{ - bra0.n_shells(), bra1.n_shells(), ket0.n_shells(), ket1.n_shells()}; + std::array n_shells{bra0.n_shells(), bra1.n_shells(), + ket0.n_shells(), ket1.n_shells()}; std::array shell_i{0, 0, 0, 0}; std::array prim_offset{0, 0, 0, 0}; std::array ao_offsets{0, 0, 0, 0}; @@ -176,35 +174,35 @@ MODULE_RUN(PrimitiveErrorModel) { // We'll collect the shells into this array std::array shells; - for(shell_i[0] = 0 ; shell_i[0] < n_shells[0]; ++shell_i[0]){ + for(shell_i[0] = 0; shell_i[0] < n_shells[0]; ++shell_i[0]) { shells[0] = bra0.shell(shell_i[0]); - + prim_offset[1] = 0; - ao_offsets[1] = 0; - for(shell_i[1] = 0; shell_i[1] < n_shells[1]; ++shell_i[1]){ + ao_offsets[1] = 0; + for(shell_i[1] = 0; shell_i[1] < n_shells[1]; ++shell_i[1]) { shells[1] = bra1.shell(shell_i[1]); - + prim_offset[2] = 0; - ao_offsets[2] = 0; - for(shell_i[2] = 0; shell_i[2] < n_shells[2]; ++shell_i[2]){ + ao_offsets[2] = 0; + for(shell_i[2] = 0; shell_i[2] < n_shells[2]; ++shell_i[2]) { shells[2] = ket0.shell(shell_i[2]); - + prim_offset[3] = 0; - ao_offsets[3] = 0; - for(shell_i[3] = 0; shell_i[3] < n_shells[3]; ++shell_i[3]){ + ao_offsets[3] = 0; + for(shell_i[3] = 0; shell_i[3] < n_shells[3]; ++shell_i[3]) { shells[3] = ket1.shell(shell_i[3]); - auto shell_error = - shell_block_error(shells, prim_offset, Q_ab, Q_cd, K_ab, K_cd, tol); - fill_ao_block(shells, ao_offsets, buffer, shell_error); - + auto shell_error = shell_block_error( + shells, prim_offset, Q_ab, Q_cd, K_ab, K_cd, tol); + fill_ao_block(shells, ao_offsets, buffer, shell_error); + prim_offset[3] += shells[3].n_primitives(); ao_offsets[3] += shells[3].size(); } // End loop over shell_i[3] prim_offset[2] += shells[2].n_primitives(); ao_offsets[2] += shells[2].size(); - }// End loop over shell_i[2] + } // End loop over shell_i[2] prim_offset[1] += shells[1].n_primitives(); ao_offsets[1] += shells[1].size(); @@ -214,10 +212,9 @@ MODULE_RUN(PrimitiveErrorModel) { ao_offsets[0] += shells[0].size(); } // End loop over shell_i[0] - simde::type::tensor error(shape, std::move(buffer)); auto result = results(); return pt::wrap_results(result, error); } -} \ No newline at end of file +} // namespace integrals::libint \ No newline at end of file diff --git a/tests/cxx/unit/integrals/ao_integrals/ao_integrals_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/ao_integrals_driver.cpp index d3ac55f2..2003d808 100644 --- a/tests/cxx/unit/integrals/ao_integrals/ao_integrals_driver.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/ao_integrals_driver.cpp @@ -14,9 +14,10 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" using simde::type::tensor; +using namespace integrals::testing; namespace { @@ -51,22 +52,20 @@ TEST_CASE("AOIntegralsDriver") { chemist::braket::BraKet; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("AO integral driver")); auto& mod = mm.at("AO integral driver"); // Get basis set - auto mol = test::h2_molecule(); - auto aobs = test::h2_sto3g_basis_set(); + auto mol = h2_molecule(); + auto aobs = h2_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); // Operator Inputs simde::type::electron e; - auto rho = test::h2_density(); + auto rho = h2_density(); SECTION("Calling Kinetic") { auto& tmod = mm.at("Kinetic"); diff --git a/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp b/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp index b49e4c48..031c59ff 100644 --- a/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp @@ -14,19 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Coulomb Metric") { using test_pt = simde::ERI2; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Coulomb Metric")); // Get basis set - auto mol = test::h2_molecule(); - auto aobs = test::h2_sto3g_basis_set(); + auto mol = h2_molecule(); + auto aobs = h2_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); @@ -40,7 +40,7 @@ TEST_CASE("Coulomb Metric") { // Call module auto T = mm.at("Coulomb Metric").run_as(braket); - auto t = test::eigen_tensor<2>(T.buffer()); + auto t = eigen_tensor<2>(T.buffer()); REQUIRE(t(0, 0) == Catch::Approx(0.15824726).margin(1E-6)); REQUIRE(t(0, 1) == Catch::Approx(0.0).margin(1E-6)); REQUIRE(t(1, 0) == Catch::Approx(-0.23097095).margin(1E-6)); diff --git a/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp b/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp index f58bbd21..7a36542e 100644 --- a/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp @@ -14,19 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Density Fitting Integral") { using test_pt = simde::ERI3; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Density Fitting Integral")); // Get basis set - auto mol = test::h2_molecule(); - auto aobs = test::h2_sto3g_basis_set(); + auto mol = h2_molecule(); + auto aobs = h2_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); @@ -41,7 +41,7 @@ TEST_CASE("Density Fitting Integral") { // Call module auto T = mm.at("Density Fitting Integral").run_as(braket); - auto t = test::eigen_tensor<3>(T.buffer()); + auto t = eigen_tensor<3>(T.buffer()); REQUIRE(t(0, 0, 0) == Catch::Approx(0.81362039).margin(1E-6)); REQUIRE(t(0, 0, 1) == Catch::Approx(0.31266336).margin(1E-6)); REQUIRE(t(0, 1, 0) == Catch::Approx(0.31266336).margin(1E-6)); diff --git a/tests/cxx/unit/integrals/ao_integrals/j_density_fitted.cpp b/tests/cxx/unit/integrals/ao_integrals/j_density_fitted.cpp index 17938991..ae01db8c 100644 --- a/tests/cxx/unit/integrals/ao_integrals/j_density_fitted.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/j_density_fitted.cpp @@ -14,25 +14,25 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Density Fitted J builder") { using pt = simde::aos_j_e_aos; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Density Fitted J builder")); // Get basis set - auto mol = test::h2_molecule(); - auto aobs = test::h2_sto3g_basis_set(); + auto mol = h2_molecule(); + auto aobs = h2_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); // Make Operator - simde::type::j_e_type op(simde::type::electron{}, test::h2_density()); + simde::type::j_e_type op(simde::type::electron{}, h2_density()); // Make BraKet Input chemist::braket::BraKet braket(aos, op, aos); @@ -41,7 +41,7 @@ TEST_CASE("Density Fitted J builder") { mm.change_input("Density Fitted J builder", "Auxiliary Basis Set", aos); const auto& T = mm.at("Density Fitted J builder").run_as(braket); - auto t = test::eigen_tensor<2>(T.buffer()); + auto t = eigen_tensor<2>(T.buffer()); REQUIRE(t(0, 0) == Catch::Approx(0.50515668).margin(1E-6)); REQUIRE(t(0, 1) == Catch::Approx(0.22344536).margin(1E-6)); REQUIRE(t(1, 0) == Catch::Approx(0.22344536).margin(1E-6)); diff --git a/tests/cxx/unit/integrals/ao_integrals/j_four_center.cpp b/tests/cxx/unit/integrals/ao_integrals/j_four_center.cpp index fe882f9b..ff13bac6 100644 --- a/tests/cxx/unit/integrals/ao_integrals/j_four_center.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/j_four_center.cpp @@ -14,25 +14,25 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Four center J builder") { using pt = simde::aos_j_e_aos; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Four center J builder")); // Get basis set - auto mol = test::h2_molecule(); - auto aobs = test::h2_sto3g_basis_set(); + auto mol = h2_molecule(); + auto aobs = h2_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); // Make Operator - simde::type::j_e_type op(simde::type::electron{}, test::h2_density()); + simde::type::j_e_type op(simde::type::electron{}, h2_density()); // Make BraKet Input chemist::braket::BraKet braket(aos, op, aos); @@ -40,7 +40,7 @@ TEST_CASE("Four center J builder") { // Call module const auto& T = mm.at("Four center J builder").run_as(braket); - auto t = test::eigen_tensor<2>(T.buffer()); + auto t = eigen_tensor<2>(T.buffer()); REQUIRE(t(0, 0) == Catch::Approx(0.56044143).margin(1E-6)); REQUIRE(t(0, 1) == Catch::Approx(0.24704427).margin(1E-6)); REQUIRE(t(1, 0) == Catch::Approx(0.24704427).margin(1E-6)); diff --git a/tests/cxx/unit/integrals/ao_integrals/k_density_fitted.cpp b/tests/cxx/unit/integrals/ao_integrals/k_density_fitted.cpp index b08664c7..fe7a8e1a 100644 --- a/tests/cxx/unit/integrals/ao_integrals/k_density_fitted.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/k_density_fitted.cpp @@ -14,25 +14,25 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Density Fitted K builder") { using pt = simde::aos_k_e_aos; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Density Fitted K builder")); // Get basis set - auto mol = test::h2_molecule(); - auto aobs = test::h2_sto3g_basis_set(); + auto mol = h2_molecule(); + auto aobs = h2_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); // Make Operator - simde::type::k_e_type op(simde::type::electron{}, test::h2_density()); + simde::type::k_e_type op(simde::type::electron{}, h2_density()); // Make BraKet Input chemist::braket::BraKet braket(aos, op, aos); @@ -41,7 +41,7 @@ TEST_CASE("Density Fitted K builder") { mm.change_input("Density Fitted K builder", "Auxiliary Basis Set", aos); const auto& T = mm.at("Density Fitted K builder").run_as(braket); - auto t = test::eigen_tensor<2>(T.buffer()); + auto t = eigen_tensor<2>(T.buffer()); REQUIRE(t(0, 0) == Catch::Approx(0.40594955).margin(1E-6)); REQUIRE(t(0, 1) == Catch::Approx(0.32265250).margin(1E-6)); REQUIRE(t(1, 0) == Catch::Approx(0.32265250).margin(1E-6)); diff --git a/tests/cxx/unit/integrals/ao_integrals/k_four_center.cpp b/tests/cxx/unit/integrals/ao_integrals/k_four_center.cpp index d5010ab2..85f06dbc 100644 --- a/tests/cxx/unit/integrals/ao_integrals/k_four_center.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/k_four_center.cpp @@ -14,25 +14,25 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Four center K builder") { using pt = simde::aos_k_e_aos; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Four center K builder")); // Get basis set - auto mol = test::h2_molecule(); - auto aobs = test::h2_sto3g_basis_set(); + auto mol = h2_molecule(); + auto aobs = h2_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); // Make Operator - simde::type::k_e_type op(simde::type::electron{}, test::h2_density()); + simde::type::k_e_type op(simde::type::electron{}, h2_density()); // Make BraKet Input chemist::braket::BraKet braket(aos, op, aos); @@ -40,7 +40,7 @@ TEST_CASE("Four center K builder") { // Call module const auto& T = mm.at("Four center K builder").run_as(braket); - auto t = test::eigen_tensor<2>(T.buffer()); + auto t = eigen_tensor<2>(T.buffer()); REQUIRE(t(0, 0) == Catch::Approx(0.45617623).margin(1E-6)); REQUIRE(t(0, 1) == Catch::Approx(0.35130947).margin(1E-6)); REQUIRE(t(1, 0) == Catch::Approx(0.35130947).margin(1E-6)); diff --git a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp index e7a9b722..398d6cec 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp @@ -14,7 +14,9 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; using namespace tensorwrapper; @@ -59,8 +61,8 @@ TEST_CASE("UQ Driver") { mm.change_input("ERI4", "Threshold", 1.0e-6); // Get basis set - auto mol = test::h2_molecule(); - auto aobs = test::h2_sto3g_basis_set(); + auto mol = h2_molecule(); + auto aobs = h2_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); diff --git a/tests/cxx/unit/integrals/libint/analytic_error.cpp b/tests/cxx/unit/integrals/libint/analytic_error.cpp index 1c95825f..1caaa132 100644 --- a/tests/cxx/unit/integrals/libint/analytic_error.cpp +++ b/tests/cxx/unit/integrals/libint/analytic_error.cpp @@ -1,16 +1,26 @@ -#include "test_error.hpp" +#include "../testing/testing.hpp" #include using eri4_pt = simde::ERI4; using pt = integrals::property_types::Uncertainty; -using namespace integrals::libint::test; +using namespace integrals::testing; +using tensorwrapper::operations::approximately_equal; +namespace { -TEST_CASE("AnalyticError") { - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); +template +auto corr_error() { + FloatType diff = -0.0000000054867100; + tensorwrapper::shape::Smooth shape({1, 1, 1, 1}); + std::vector buffer(shape.size(), diff); + tensorwrapper::buffer::Contiguous cont(std::move(buffer), shape); + return simde::type::tensor(std::move(cont), shape); +} + +} // namespace +TEST_CASE("AnalyticError") { + auto mm = initialize_integrals(); auto& mod = mm.at("Analytic Error"); simde::type::v_ee_type v_ee{}; @@ -24,11 +34,7 @@ TEST_CASE("AnalyticError") { chemist::braket::BraKet mnls(bra, v_ee, ket); float_type tol = 1.0e-10; auto error = mod.run_as(mnls, tol); - tensorwrapper::shape::Smooth shape({1, 1, 1, 1}); - auto benchmark = test::make_tw_buffer(0.0002263495626484, shape); - auto screened = test::make_tw_buffer(0.0002263440759384, shape); - tensorwrapper::buffer::Contiguous corr(benchmark); - corr("m,n,l,s") = screened("m,n,l,s") - benchmark("m,n,l,s"); - REQUIRE(error.buffer().approximately_equal(corr, 1.0e-12)); + auto corr = corr_error(); + REQUIRE(approximately_equal(error, corr, 1.0e-12)); } } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp index 65df1b71..52134ed5 100644 --- a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp +++ b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp @@ -13,50 +13,55 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#include "testing.hpp" +#include "../testing/testing.hpp" #include +using namespace integrals::testing; +using tensorwrapper::operations::approximately_equal; + using pt = integrals::property_types::PrimitivePairEstimator; namespace { -// std::vector h2_h2_corr{ -// 0.023817430147884473, 0.08261663936778646, 0.06861998972363427, -// 1.4555047643750448e-05, 0.00844595180076824, 0.03423471319097009, -// 0.08261663936778646, 0.28657621993836907, 0.23802538347833696, -// 0.00844595180076824, 0.07444365105885908, 0.13404298325096972, -// 0.06861998972363427, 0.23802538347833696, 0.19769987611740356, -// 0.03423471319097009, 0.13404298325096972, 0.1372685140907819, -// 1.4555047643750448e-05, 0.00844595180076824, 0.03423471319097009, -// 0.023817430147884473, 0.08261663936778646, 0.06861998972363427, -// 0.00844595180076824, 0.07444365105885908, 0.13404298325096972, -// 0.08261663936778646, 0.28657621993836907, 0.23802538347833696, -// 0.03423471319097009, 0.13404298325096972, 0.1372685140907819, -// 0.06861998972363427, 0.23802538347833696, 0.19769987611740356}; - -// } // namespace +template +auto corr_k01() { + std::vector buffer{3.693e-38, 4.76823e-13, 8.9425e-06, + 4.76823e-13, 1.73977e-08, 3.39937e-05, + 8.9425e-06, 3.39937e-05, 0.000112915}; + tensorwrapper::shape::Smooth shape({3, 3}); + tensorwrapper::buffer::Contiguous cont(std::move(buffer), shape); + return simde::type::tensor(std::move(cont), shape); +} + +template +auto corr_k23() { + std::vector buffer{4.07419e-12, 5.0571e-05, 0.00250317, + 5.0571e-05, 0.000964263, 0.00356585, + 0.00250317, 0.00356585, 0.00217062}; + + tensorwrapper::shape::Smooth shape({3, 3}); + tensorwrapper::buffer::Contiguous cont(std::move(buffer), shape); + return simde::type::tensor(std::move(cont), shape); +} + +} // namespace TEST_CASE("BlackBoxPrimitiveEstimator") { - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - auto& mod = mm.at("Black Box Primitive Pair Estimator"); - auto h2_basis = test::h2_sto3g_basis_set(); - - SECTION("H2 STO-3G") { - const auto& Kij = mod.run_as(h2_basis, h2_basis); - tensorwrapper::shape::Smooth shape{6, 6}; - tensorwrapper::buffer::Contiguous corr(h2_h2_corr, shape); - REQUIRE(Kij.buffer().approximately_equal(corr, 1e-8)); + auto mm = initialize_integrals(); + auto& mod = mm.at("Black Box Primitive Pair Estimator"); + auto&& [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + + using float_type = double; + + SECTION("K(bra0, bra1)") { + auto K01 = mod.run_as(bra0, bra1); + auto corr = corr_k01(); + REQUIRE(approximately_equal(K01, corr, 1E-8)); } - // Tests bra != ket - SECTION("Block of H2 STO-3G") { - simde::type::ao_basis_set h0_basis; - h0_basis.add_center(h2_basis.at(0)); - const auto& Kij = mod.run_as(h0_basis, h2_basis); - tensorwrapper::shape::Smooth shape{3, 6}; - std::vector buffer(h2_h2_corr.begin(), h2_h2_corr.begin() + 18); - tensorwrapper::buffer::Contiguous corr(buffer, shape); - REQUIRE(Kij.buffer().approximately_equal(corr, 1e-8)); + SECTION("K(ket0, ket1)") { + auto K23 = mod.run_as(ket0, ket1); + auto corr = corr_k23(); + REQUIRE(approximately_equal(K23, corr, 1E-8)); } } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp index 5a6c70e9..252b1e4c 100644 --- a/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp +++ b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -1,31 +1,29 @@ -/* - * 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 "test_error.hpp" -#include +// /* +// * 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 "../testing/testing.hpp" +// #include -using pt = integrals::property_types::PrimitivePairEstimator; -using namespace integrals::libint::test; -TEST_CASE("CauchySchwarzPrimitiveEstimator") { - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); - auto& mod = mm.at("CauchySchwarz Estimator"); - auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); +// using pt = integrals::property_types::PrimitivePairEstimator; +// using namespace integrals::testing; +// TEST_CASE("CauchySchwarzPrimitiveEstimator") { +// auto mm = initialize_integrals(); +// auto& mod = mm.at("CauchySchwarz Estimator"); +// auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); - auto Q_ab = mod.run_as(bra0, bra1); - auto Q_cd = mod.run_as(ket0, ket1); - // std::cout << Q_ab << std::endl; -} \ No newline at end of file +// auto Q_ab = mod.run_as(bra0, bra1); +// auto Q_cd = mod.run_as(ket0, ket1); +// // std::cout << Q_ab << std::endl; +// } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/detail_/test_make_libint_basis_set.cpp b/tests/cxx/unit/integrals/libint/detail_/test_make_libint_basis_set.cpp index 991b084e..1e4d2b57 100644 --- a/tests/cxx/unit/integrals/libint/detail_/test_make_libint_basis_set.cpp +++ b/tests/cxx/unit/integrals/libint/detail_/test_make_libint_basis_set.cpp @@ -18,11 +18,13 @@ #include "libint_basis_set_water.hpp" #undef DEPRECATED // Must be last due to conflicting macros -#include "../../testing.hpp" +#include "../../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("make_libint_basis_set") { using integrals::libint::detail_::make_libint_basis_set; - auto aobs = test::water_sto3g_basis_set(); + auto aobs = water_sto3g_basis_set(); auto libint_bs = make_libint_basis_set(aobs); auto libint_corr = test::water_basis_set(); REQUIRE(libint_bs == libint_corr); diff --git a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp index 42a45938..4df662db 100644 --- a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -1,52 +1,50 @@ -#include "test_error.hpp" -#include - -using eri4_pt = simde::ERI4; -using pt = integrals::property_types::Uncertainty; - -using namespace integrals::libint::test; - -TEST_CASE("PrimitiveErrorModel") { - pluginplay::ModuleManager mm; - integrals::load_modules(mm); - integrals::set_defaults(mm); - - auto& mod = mm.at("Primitive Error Model"); - - simde::type::v_ee_type v_ee{}; - using float_type = double; - SECTION("H2") { - auto h2_aos = test::h2_sto3g_basis_set(); - simde::type::aos_squared h2_aos2(h2_aos, h2_aos); - chemist::braket::BraKet mnls(h2_aos2, v_ee, h2_aos2); - auto error = mod.run_as(mnls, 1.0e-10); - - tensorwrapper::shape::Smooth shape({2, 2, 2, 2}); - std::vector data(shape.size()); - tensorwrapper::buffer::Contiguous buffer(std::move(data), shape); - REQUIRE(buffer.approximately_equal(error.buffer(), 1.0e-12)); - } - - SECTION("H2 2") { - float_type tol = 1.0e-10; - auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); - auto value = manual_contract_shell(bra0, bra1, ket0, ket1, mm); - std::cout << "Screened value: " << value << std::endl; - - simde::type::aos_squared bra(bra0, bra1); - simde::type::aos_squared ket(ket0, ket1); - chemist::braket::BraKet mnls(bra, v_ee, ket); - - auto error = mod.run_as(mnls, tol); - auto corr = mm.at("Analytic Error").run_as(mnls, tol); - - auto eri4_mod = mm.at("ERI4").unlocked_copy(); - eri4_mod.change_input("Threshold", float_type(1E-16)); - auto i_value = eri4_mod.run_as(mnls); - std::cout << "Corr Screened Value: " << i_value << std::endl; - - // Our value: 0.0002263495591894 - // Libint's value: 0.0002263440759384 - // No screening: 0.0002263495626484 - } -} \ No newline at end of file +// #include "../testing/testing.hpp" +// #include + +// using eri4_pt = simde::ERI4; +// using pt = integrals::property_types::Uncertainty; + +// using namespace integrals::testing; + +// TEST_CASE("PrimitiveErrorModel") { +// auto mm = initialize_integrals(); + +// auto& mod = mm.at("Primitive Error Model"); + +// simde::type::v_ee_type v_ee{}; +// using float_type = double; +// SECTION("H2") { +// auto h2_aos = test::h2_sto3g_basis_set(); +// simde::type::aos_squared h2_aos2(h2_aos, h2_aos); +// chemist::braket::BraKet mnls(h2_aos2, v_ee, h2_aos2); +// auto error = mod.run_as(mnls, 1.0e-10); + +// tensorwrapper::shape::Smooth shape({2, 2, 2, 2}); +// std::vector data(shape.size()); +// tensorwrapper::buffer::Contiguous buffer(std::move(data), shape); +// REQUIRE(buffer.approximately_equal(error.buffer(), 1.0e-12)); +// } + +// SECTION("H2 2") { +// float_type tol = 1.0e-10; +// auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); +// auto value = manual_contract_shell(bra0, bra1, ket0, ket1, mm); +// std::cout << "Screened value: " << value << std::endl; + +// simde::type::aos_squared bra(bra0, bra1); +// simde::type::aos_squared ket(ket0, ket1); +// chemist::braket::BraKet mnls(bra, v_ee, ket); + +// auto error = mod.run_as(mnls, tol); +// auto corr = mm.at("Analytic Error").run_as(mnls, tol); + +// auto eri4_mod = mm.at("ERI4").unlocked_copy(); +// eri4_mod.change_input("Threshold", float_type(1E-16)); +// auto i_value = eri4_mod.run_as(mnls); +// std::cout << "Corr Screened Value: " << i_value << std::endl; + +// // Our value: 0.0002263495591894 +// // Libint's value: 0.0002263440759384 +// // No screening: 0.0002263495626484 +// } +// } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/test_arbitrary_operator.cpp b/tests/cxx/unit/integrals/libint/test_arbitrary_operator.cpp index 9c1e5680..be7667f9 100644 --- a/tests/cxx/unit/integrals/libint/test_arbitrary_operator.cpp +++ b/tests/cxx/unit/integrals/libint/test_arbitrary_operator.cpp @@ -14,9 +14,11 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" #include "integrals/uncertain_types.hpp" +using namespace integrals::testing; + using udouble = integrals::type::uncertain_double; constexpr bool has_sigma = integrals::type::has_sigma(); @@ -44,15 +46,14 @@ TEST_CASE("OperatorBase") { using op_t = simde::type::v_ee_type; using op_base_t = simde::type::op_base_type; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Evaluate 2-Index BraKet")); REQUIRE(mm.count("Evaluate 3-Index BraKet")); REQUIRE(mm.count("Evaluate 4-Index BraKet")); // Get basis set - auto mol = test::water_molecule(); - auto aobs = test::water_sto3g_basis_set(); + auto mol = water_molecule(); + auto aobs = water_sto3g_basis_set(); // Make AOS object aos_t aos(aobs); @@ -76,9 +77,9 @@ TEST_CASE("OperatorBase") { auto T = mod.run_as(braket); // Check output - REQUIRE(test::trace<2>(T.buffer()) == + REQUIRE(trace<2>(T.buffer()) == Catch::Approx(124.7011973877891364).margin(1.0e-16)); - REQUIRE(test::norm<2>(T.buffer()) == + REQUIRE(norm<2>(T.buffer()) == Catch::Approx(90.2562579028763707).margin(1.0e-16)); } @@ -88,13 +89,13 @@ TEST_CASE("OperatorBase") { auto T = mod.run_as(braket); // Check output - REQUIRE(unwrap_mean(test::trace<2, udouble>(T.buffer())) == + REQUIRE(unwrap_mean(trace<2, udouble>(T.buffer())) == Catch::Approx(124.7011973877891364).margin(1.0e-16)); - REQUIRE(unwrap_sd(test::trace<2, udouble>(T.buffer())) == + REQUIRE(unwrap_sd(trace<2, udouble>(T.buffer())) == Catch::Approx(7e-16).margin(1.0e-16)); - REQUIRE(unwrap_mean(test::norm<2, udouble>(T.buffer())) == + REQUIRE(unwrap_mean(norm<2, udouble>(T.buffer())) == Catch::Approx(90.2562579028763707).margin(1.0e-16)); - REQUIRE(unwrap_sd(test::norm<2, udouble>(T.buffer())) == + REQUIRE(unwrap_sd(norm<2, udouble>(T.buffer())) == Catch::Approx(3e-16).margin(1.0e-16)); } } @@ -113,9 +114,9 @@ TEST_CASE("OperatorBase") { auto T = mod.run_as(braket); // Check output - REQUIRE(test::trace<3>(T.buffer()) == + REQUIRE(trace<3>(T.buffer()) == Catch::Approx(16.8245948391706577).margin(1.0e-16)); - REQUIRE(test::norm<3>(T.buffer()) == + REQUIRE(norm<3>(T.buffer()) == Catch::Approx(20.6560572032543597).margin(1.0e-16)); } @@ -127,13 +128,13 @@ TEST_CASE("OperatorBase") { // Check output auto& t = T.buffer(); - REQUIRE(unwrap_mean(test::trace<3, udouble>(t)) == + REQUIRE(unwrap_mean(trace<3, udouble>(t)) == Catch::Approx(16.8245948391706577).margin(1.0e-16)); - REQUIRE(unwrap_sd(test::trace<3, udouble>(t)) == + REQUIRE(unwrap_sd(trace<3, udouble>(t)) == Catch::Approx(7e-16).margin(1.0e-16)); - REQUIRE(unwrap_mean(test::norm<3, udouble>(t)) == + REQUIRE(unwrap_mean(norm<3, udouble>(t)) == Catch::Approx(20.6560572032543597).margin(1.0e-16)); - REQUIRE(unwrap_sd(test::norm<3, udouble>(t)) == + REQUIRE(unwrap_sd(norm<3, udouble>(t)) == Catch::Approx(7e-16).margin(1.0e-16)); } } @@ -155,9 +156,9 @@ TEST_CASE("OperatorBase") { // Check output auto& t = T.buffer(); - REQUIRE(test::trace<4>(t) == + REQUIRE(trace<4>(t) == Catch::Approx(9.7919608941952063).margin(1.0e-16)); - REQUIRE(test::norm<4>(t) == + REQUIRE(norm<4>(t) == Catch::Approx(7.7796143419802553).margin(1.0e-16)); } @@ -169,13 +170,13 @@ TEST_CASE("OperatorBase") { // Check output auto& t = T.buffer(); - REQUIRE(unwrap_mean(test::trace<4, udouble>(t)) == + REQUIRE(unwrap_mean(trace<4, udouble>(t)) == Catch::Approx(9.7919608941952063).margin(1.0e-16)); - REQUIRE(unwrap_sd(test::trace<4, udouble>(t)) == + REQUIRE(unwrap_sd(trace<4, udouble>(t)) == Catch::Approx(7e-16).margin(1.0e-16)); - REQUIRE(unwrap_mean(test::norm<4, udouble>(t)) == + REQUIRE(unwrap_mean(norm<4, udouble>(t)) == Catch::Approx(7.7796143419802553).margin(1.0e-16)); - REQUIRE(unwrap_sd(test::norm<4, udouble>(t)) == + REQUIRE(unwrap_sd(norm<4, udouble>(t)) == Catch::Approx(11e-16).margin(1.0e-16)); } } diff --git a/tests/cxx/unit/integrals/libint/test_eri2.cpp b/tests/cxx/unit/integrals/libint/test_eri2.cpp index b98afd13..3e58f870 100644 --- a/tests/cxx/unit/integrals/libint/test_eri2.cpp +++ b/tests/cxx/unit/integrals/libint/test_eri2.cpp @@ -14,18 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("ERI2") { using test_pt = simde::ERI2; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("ERI2")); // Get basis set - auto mol = test::water_molecule(); - auto aobs = test::water_sto3g_basis_set(); + auto mol = water_molecule(); + auto aobs = water_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); @@ -41,8 +42,6 @@ TEST_CASE("ERI2") { // Check output auto& t = T.buffer(); - REQUIRE(test::trace<2>(t) == - Catch::Approx(124.7011973877891364).margin(1.0e-16)); - REQUIRE(test::norm<2>(t) == - Catch::Approx(90.2562579028763707).margin(1.0e-16)); + REQUIRE(trace<2>(t) == Catch::Approx(124.7011973877891364).margin(1.0e-16)); + REQUIRE(norm<2>(t) == Catch::Approx(90.2562579028763707).margin(1.0e-16)); } diff --git a/tests/cxx/unit/integrals/libint/test_eri3.cpp b/tests/cxx/unit/integrals/libint/test_eri3.cpp index 43d4287c..cc00a10b 100644 --- a/tests/cxx/unit/integrals/libint/test_eri3.cpp +++ b/tests/cxx/unit/integrals/libint/test_eri3.cpp @@ -14,18 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("ERI3") { using test_pt = simde::ERI3; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("ERI3")); // Get basis set - auto mol = test::water_molecule(); - auto aobs = test::water_sto3g_basis_set(); + auto mol = water_molecule(); + auto aobs = water_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); @@ -42,8 +43,6 @@ TEST_CASE("ERI3") { // Check output auto& t = T.buffer(); - REQUIRE(test::trace<3>(t) == - Catch::Approx(16.8245948391706577).margin(1.0e-16)); - REQUIRE(test::norm<3>(t) == - Catch::Approx(20.6560572032543597).margin(1.0e-16)); + REQUIRE(trace<3>(t) == Catch::Approx(16.8245948391706577).margin(1.0e-16)); + REQUIRE(norm<3>(t) == Catch::Approx(20.6560572032543597).margin(1.0e-16)); } diff --git a/tests/cxx/unit/integrals/libint/test_eri4.cpp b/tests/cxx/unit/integrals/libint/test_eri4.cpp index 43881954..16957d5a 100644 --- a/tests/cxx/unit/integrals/libint/test_eri4.cpp +++ b/tests/cxx/unit/integrals/libint/test_eri4.cpp @@ -14,18 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("ERI4") { using test_pt = simde::ERI4; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("ERI4")); // Get basis set - auto mol = test::water_molecule(); - auto aobs = test::water_sto3g_basis_set(); + auto mol = water_molecule(); + auto aobs = water_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); @@ -42,8 +43,6 @@ TEST_CASE("ERI4") { // Check output auto& t = T.buffer(); - REQUIRE(test::trace<4>(t) == - Catch::Approx(9.7919608941952063).margin(1.0e-16)); - REQUIRE(test::norm<4>(t) == - Catch::Approx(7.7796143419802553).margin(1.0e-16)); + REQUIRE(trace<4>(t) == Catch::Approx(9.7919608941952063).margin(1.0e-16)); + REQUIRE(norm<4>(t) == Catch::Approx(7.7796143419802553).margin(1.0e-16)); } diff --git a/tests/cxx/unit/integrals/libint/test_error.hpp b/tests/cxx/unit/integrals/libint/test_error.hpp deleted file mode 100644 index 32cfa4e1..00000000 --- a/tests/cxx/unit/integrals/libint/test_error.hpp +++ /dev/null @@ -1,123 +0,0 @@ -#pragma once -#include "testing.hpp" - -namespace integrals::libint::test { - -/* After playing around, I found that if you set up a H2 dimer so that the - H2 molecules are 3.0 Angstroms apart, then the (03|12) integral has a - relatively large error (about 5.5E-9 using a tolerance of 1.0E-10). This - function prepares 4, 1-shell AO basis sets so that you can easily evaluate - that element. - */ -inline auto get_h2_dimer_0312_bases() { - simde::type::ao_basis_set bra0, bra1, ket0, ket1; - auto h2_2_aos = ::test::h2_sto3g_basis_set(); - - // Centers 0 and 1 stay put, center 2 is 0 translated, and 3 is 1 translated - bra0.add_center(h2_2_aos[0]); // 0 - bra1.add_center(h2_2_aos[1]); // 3 - ket0.add_center(h2_2_aos[1]); // 1 - ket1.add_center(h2_2_aos[0]); // 2 - - const double r = 3.0; - auto r0 = h2_2_aos[0].center(); - auto r1 = h2_2_aos[1].center(); - - bra1[0].center().x() = r1.x() + r; - bra1[0].center().y() = r1.y() + r; - bra1[0].center().z() = r1.z() + r; - ket1[0].center().x() = r0.x() + r; - ket1[0].center().y() = r0.y() + r; - ket1[0].center().z() = r0.z() + r; - return std::make_tuple(bra0, bra1, ket0, ket1); -} - -template -double black_box_pair_metric(const T& prim_i, const T& prim_j) { - auto ci = prim_i.coefficient(); - auto cj = prim_j.coefficient(); - auto zeta_i = prim_i.exponent(); - auto zeta_j = prim_j.exponent(); - auto dx = prim_i.center().x() - prim_j.center().x(); - auto dy = prim_i.center().y() - prim_j.center().y(); - auto dz = prim_i.center().z() - prim_j.center().z(); - auto dr = std::sqrt(dx * dx + dy * dy + dz * dz); - auto exponent = (-zeta_i * zeta_j / (zeta_i + zeta_j)) * dr * dr; - return ci * cj * std::exp(exponent); -} - -/* Given four 1 shell basis sets, this function will create the shell quartet - in the AO basis set by contracting the primitive integrals with the - primitive coefficients. This is a sanity check more than anything else. - */ -template -inline auto manual_contract_shell(const BasisType& bra0, const BasisType& bra1, - const BasisType& ket0, const BasisType& ket1, - pluginplay::ModuleManager& mm) { - auto dec_mod = mm.at("Decontract Basis Set"); - using dec_pt = integrals::property_types::DecontractBasisSet; - auto bra0_prims = dec_mod.run_as(bra0); - auto bra1_prims = dec_mod.run_as(bra1); - auto ket0_prims = dec_mod.run_as(ket0); - auto ket1_prims = dec_mod.run_as(ket1); - simde::type::v_ee_type v_ee{}; - simde::type::aos_squared bra_prims(bra0_prims, bra1_prims); - simde::type::aos_squared ket_prims(ket0_prims, ket1_prims); - chemist::braket::BraKet mnls(bra_prims, v_ee, ket_prims); - auto ints = mm.at("ERI4").run_as(mnls); - // auto K_mod = mm.at("Black Box Primitive Pair Estimator"); - // using prim_pt = integrals::property_types::PrimitivePairEstimator; - // auto K_ab = K_mod.run_as(bra0, bra1); - // auto K_cd = K_mod.run_as(ket0, ket1); - - const auto& buffer = tensorwrapper::buffer::make_contiguous(ints.buffer()); - // const auto& K01 = - // tensorwrapper::buffer::make_contiguous(K_ab.buffer()); const auto& K23 = - // tensorwrapper::buffer::make_contiguous(K_cd.buffer()); - - // TODO: expand to higher angular momentum - double value = 0.0; - double tol = 1e-10; - double error = 0.0; - using wtf::fp::float_cast; - for(std::size_t i = 0; i < bra0_prims.n_primitives(); ++i) { - const auto prim_i = bra0.primitive(i); - auto ci = prim_i.coefficient(); - - for(std::size_t j = 0; j < bra1_prims.n_primitives(); ++j) { - const auto prim_j = bra1.primitive(j); - auto cj = prim_j.coefficient(); - - auto kij = black_box_pair_metric(prim_i, prim_j); - - bool is01good = kij > tol; - - for(std::size_t k = 0; k < ket0_prims.n_primitives(); ++k) { - const auto prim_k = ket0.primitive(k); - auto ck = prim_k.coefficient(); - - for(std::size_t l = 0; l < ket1_prims.n_primitives(); ++l) { - const auto prim_l = ket1.primitive(l); - auto cl = prim_l.coefficient(); - - auto kkl = black_box_pair_metric(prim_k, prim_l); - bool is23good = kkl > tol; - bool both_good = (kij * kkl) > tol; - - auto erased = buffer.get_elem({i, j, k, l}); - auto val = wtf::fp::float_cast(erased); - double Iijk = ci * cj * ck * cl * val; - if(is01good && is23good && both_good) { - value += Iijk; - } else { - error += std::fabs(Iijk); - } - } - } - } - } - std::cout << "Computed error: " << simde::type::tensor(error) << std::endl; - return simde::type::tensor(value); -} - -} // namespace integrals::libint::test \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/test_kinetic.cpp b/tests/cxx/unit/integrals/libint/test_kinetic.cpp index 9f000c25..6eace9b7 100644 --- a/tests/cxx/unit/integrals/libint/test_kinetic.cpp +++ b/tests/cxx/unit/integrals/libint/test_kinetic.cpp @@ -14,18 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Kinetic") { using test_pt = simde::aos_t_e_aos; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Kinetic")); // Get basis set - auto mol = test::water_molecule(); - auto aobs = test::water_sto3g_basis_set(); + auto mol = water_molecule(); + auto aobs = water_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); @@ -41,8 +42,6 @@ TEST_CASE("Kinetic") { // Check output auto& t = T.buffer(); - REQUIRE(test::trace<2>(t) == - Catch::Approx(38.9175852621874441).margin(1.0e-16)); - REQUIRE(test::norm<2>(t) == - Catch::Approx(29.3665362218072552).margin(1.0e-16)); + REQUIRE(trace<2>(t) == Catch::Approx(38.9175852621874441).margin(1.0e-16)); + REQUIRE(norm<2>(t) == Catch::Approx(29.3665362218072552).margin(1.0e-16)); } diff --git a/tests/cxx/unit/integrals/libint/test_nuclear.cpp b/tests/cxx/unit/integrals/libint/test_nuclear.cpp index 9f3d8360..b792a708 100644 --- a/tests/cxx/unit/integrals/libint/test_nuclear.cpp +++ b/tests/cxx/unit/integrals/libint/test_nuclear.cpp @@ -14,18 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Nuclear") { using test_pt = simde::aos_v_en_aos; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Nuclear")); // Get basis set - auto mol = test::water_molecule(); - auto aobs = test::water_sto3g_basis_set(); + auto mol = water_molecule(); + auto aobs = water_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); @@ -41,8 +42,7 @@ TEST_CASE("Nuclear") { // Check output auto& t = T.buffer(); - REQUIRE(test::trace<2>(t) == + REQUIRE(trace<2>(t) == Catch::Approx(-111.9975421879705664).margin(1.0e-16)); - REQUIRE(test::norm<2>(t) == - Catch::Approx(66.4857539908047528).margin(1.0e-16)); + REQUIRE(norm<2>(t) == Catch::Approx(66.4857539908047528).margin(1.0e-16)); } diff --git a/tests/cxx/unit/integrals/libint/test_overlap.cpp b/tests/cxx/unit/integrals/libint/test_overlap.cpp index bece7ae3..e427e477 100644 --- a/tests/cxx/unit/integrals/libint/test_overlap.cpp +++ b/tests/cxx/unit/integrals/libint/test_overlap.cpp @@ -14,18 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" + +using namespace integrals::testing; TEST_CASE("Overlap") { using test_pt = simde::aos_s_e_aos; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Overlap")); // Get basis set - auto mol = test::water_molecule(); - auto aobs = test::water_sto3g_basis_set(); + auto mol = water_molecule(); + auto aobs = water_sto3g_basis_set(); // Make AOS object simde::type::aos aos(aobs); @@ -41,8 +42,6 @@ TEST_CASE("Overlap") { // Check output auto& t = S.buffer(); - REQUIRE(test::trace<2>(t) == - Catch::Approx(7.00000000000000266).margin(1.0e-16)); - REQUIRE(test::norm<2>(t) == - Catch::Approx(2.87134497074907324).margin(1.0e-16)); + REQUIRE(trace<2>(t) == Catch::Approx(7.00000000000000266).margin(1.0e-16)); + REQUIRE(norm<2>(t) == Catch::Approx(2.87134497074907324).margin(1.0e-16)); } diff --git a/tests/cxx/unit/integrals/testing.hpp b/tests/cxx/unit/integrals/testing/ao_bases.hpp similarity index 54% rename from tests/cxx/unit/integrals/testing.hpp rename to tests/cxx/unit/integrals/testing/ao_bases.hpp index b56a0c32..e4a619cd 100644 --- a/tests/cxx/unit/integrals/testing.hpp +++ b/tests/cxx/unit/integrals/testing/ao_bases.hpp @@ -1,80 +1,33 @@ -/* - * Copyright 2024 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 - -// For common inputs +#include "molecules.hpp" #include -// Common Catch2 includes -#include -#include - -#include - -namespace test { - -template -auto make_tw_buffer(FloatType value, tensorwrapper::shape::Smooth shape) { - std::vector data(1, value); - return tensorwrapper::buffer::Contiguous(std::move(data), shape); -} - -template -auto eigen_tensor_(const tensorwrapper::buffer::BufferBase& buffer, - std::array extents, std::index_sequence) { - using namespace tensorwrapper; - const auto pdata = buffer::get_raw_data(buffer); - using eigen_type = Eigen::Tensor; - return Eigen::TensorMap(pdata.data(), extents[Is]...); -} +namespace integrals::testing { -// Checking eigen outputs -template -auto eigen_tensor(const tensorwrapper::buffer::BufferBase& buffer) { - std::array extents; - auto shape = buffer.layout().shape().as_smooth(); - for(std::size_t i = 0; i < N; ++i) extents[i] = shape.extent(i); - return eigen_tensor_(buffer, extents, - std::make_index_sequence()); -} +inline simde::type::ao_basis_set h2_sto3g_basis_set() { + using ao_basis_t = simde::type::ao_basis_set; + using atomic_basis_t = simde::type::atomic_basis_set; + using cg_t = simde::type::contracted_gaussian; + using point_t = simde::type::point; + using doubles_t = std::vector; -template -auto trace(const tensorwrapper::buffer::BufferBase& buffer) { - auto t = eigen_tensor(buffer); - Eigen::Tensor trace = t.trace(); - return trace.coeff(); -} + auto mol = water_molecule(); + point_t r0 = mol[0].as_nucleus(); + point_t r1 = mol[1].as_nucleus(); -template -auto norm(const tensorwrapper::buffer::BufferBase& buffer) { - auto t = eigen_tensor(buffer); - Eigen::Tensor norm = t.square().sum().sqrt(); - return norm.coeff(); -} + doubles_t cs{0.1543289673, 0.5353281423, 0.4446345422}; + doubles_t es{3.425250914, 0.6239137298, 0.1688554040}; + cg_t cg0(cs.begin(), cs.end(), es.begin(), es.end(), r0); + cg_t cg1(cs.begin(), cs.end(), es.begin(), es.end(), r1); + atomic_basis_t h0("sto-3g", 1, r0); + atomic_basis_t h1("sto-3g", 1, r1); + h0.add_shell(chemist::ShellType::cartesian, 0, cg0); + h1.add_shell(chemist::ShellType::cartesian, 0, cg1); -// Inputs for Water tests -inline simde::type::molecule water_molecule() { - using atom_t = simde::type::atom; - using molecule_t = simde::type::molecule; - atom_t O{"O", 8ul, 0.0, 0.0, -0.143222342980786, 0.0}; - atom_t H1{"H", 1ul, 0.0, 1.638033502034240, 1.136556880358410, 0.0}; - atom_t H2{"H", 1ul, 0.0, -1.638033502034240, 1.136556880358410, 0.0}; - return molecule_t{O, H1, H2}; + ao_basis_t bs; + bs.add_center(h0); + bs.add_center(h1); + return bs; } inline simde::type::ao_basis_set water_sto3g_basis_set() { @@ -173,54 +126,4 @@ inline simde::type::ao_basis_set water_decontracted_sto3g_basis_set() { return bs; } -// Inputs for H2 tests -inline simde::type::molecule h2_molecule() { - using atom_t = simde::type::atom; - using molecule_t = simde::type::molecule; - atom_t H0("H", 1ul, 1836.15, 0.0, 0.0, 0.0); - atom_t H1("H", 1ul, 1836.15, 0.0, 0.0, 1.3984); - return molecule_t{H0, H1}; -} - -inline simde::type::ao_basis_set h2_sto3g_basis_set() { - using ao_basis_t = simde::type::ao_basis_set; - using atomic_basis_t = simde::type::atomic_basis_set; - using cg_t = simde::type::contracted_gaussian; - using point_t = simde::type::point; - using doubles_t = std::vector; - - // Yes this is wrong, but it's what all the tests are expecting now - auto mol = water_molecule(); - point_t r0 = mol[0].as_nucleus(); - point_t r1 = mol[1].as_nucleus(); - - doubles_t cs{0.1543289673, 0.5353281423, 0.4446345422}; - doubles_t es{3.425250914, 0.6239137298, 0.1688554040}; - cg_t cg0(cs.begin(), cs.end(), es.begin(), es.end(), r0); - cg_t cg1(cs.begin(), cs.end(), es.begin(), es.end(), r1); - atomic_basis_t h0("sto-3g", 1, r0); - atomic_basis_t h1("sto-3g", 1, r1); - h0.add_shell(chemist::ShellType::cartesian, 0, cg0); - h1.add_shell(chemist::ShellType::cartesian, 0, cg1); - - ao_basis_t bs; - bs.add_center(h0); - bs.add_center(h1); - return bs; -} - -inline auto h2_mos() { - using mos_type = simde::type::mos; - using tensor_type = typename mos_type::transform_type; - tensor_type c({{-0.565516, -1.07019}, {-0.565516, 1.07019}}); - return mos_type(simde::type::aos(h2_sto3g_basis_set()), std::move(c)); -} - -inline auto h2_density() { - using density_type = simde::type::decomposable_e_density; - typename density_type::value_type rho( - {{0.31980835, 0.31980835}, {0.31980835, 0.31980835}}); - return density_type(rho, h2_mos()); -} - -} // namespace test +} // namespace integrals::testing \ No newline at end of file diff --git a/tests/cxx/unit/integrals/testing/molecules.hpp b/tests/cxx/unit/integrals/testing/molecules.hpp new file mode 100644 index 00000000..ac87f77d --- /dev/null +++ b/tests/cxx/unit/integrals/testing/molecules.hpp @@ -0,0 +1,25 @@ +#pragma once +#include + +namespace integrals::testing { + +// Inputs for H2 tests +inline simde::type::molecule h2_molecule() { + using atom_t = simde::type::atom; + using molecule_t = simde::type::molecule; + atom_t H0("H", 1ul, 1836.15, 0.0, 0.0, 0.0); + atom_t H1("H", 1ul, 1836.15, 0.0, 0.0, 1.3984); + return molecule_t{H0, H1}; +} + +// Inputs for Water tests +inline simde::type::molecule water_molecule() { + using atom_t = simde::type::atom; + using molecule_t = simde::type::molecule; + atom_t O{"O", 8ul, 0.0, 0.0, -0.143222342980786, 0.0}; + atom_t H1{"H", 1ul, 0.0, 1.638033502034240, 1.136556880358410, 0.0}; + atom_t H2{"H", 1ul, 0.0, -1.638033502034240, 1.136556880358410, 0.0}; + return molecule_t{O, H1, H2}; +} + +} // namespace integrals::testing \ No newline at end of file diff --git a/tests/cxx/unit/integrals/testing/shell_quartets.hpp b/tests/cxx/unit/integrals/testing/shell_quartets.hpp new file mode 100644 index 00000000..2670903e --- /dev/null +++ b/tests/cxx/unit/integrals/testing/shell_quartets.hpp @@ -0,0 +1,36 @@ +#pragma once +#include "ao_bases.hpp" +#include + +namespace integrals::testing { + +/* After playing around, I found that if you set up a H2 dimer so that the + H2 molecules are 3.0 Angstroms apart, then the (03|12) integral has a + relatively large error (about 5.5E-9 using a tolerance of 1.0E-10). This + function prepares 4, 1-shell AO basis sets so that you can easily evaluate + that element. + */ +inline auto get_h2_dimer_0312_bases() { + simde::type::ao_basis_set bra0, bra1, ket0, ket1; + auto h2_2_aos = h2_sto3g_basis_set(); + + // Centers 0 and 1 stay put, center 2 is 0 translated, and 3 is 1 translated + bra0.add_center(h2_2_aos[0]); // 0 + bra1.add_center(h2_2_aos[1]); // 3 + ket0.add_center(h2_2_aos[1]); // 1 + ket1.add_center(h2_2_aos[0]); // 2 + + const double r = 3.0; + auto r0 = h2_2_aos[0].center(); + auto r1 = h2_2_aos[1].center(); + + bra1[0].center().x() = r1.x() + r; + bra1[0].center().y() = r1.y() + r; + bra1[0].center().z() = r1.z() + r; + ket1[0].center().x() = r0.x() + r; + ket1[0].center().y() = r0.y() + r; + ket1[0].center().z() = r0.z() + r; + return std::make_tuple(bra0, bra1, ket0, ket1); +} + +} // namespace integrals::testing \ No newline at end of file diff --git a/tests/cxx/unit/integrals/testing/testing.hpp b/tests/cxx/unit/integrals/testing/testing.hpp new file mode 100644 index 00000000..cb6a7e3b --- /dev/null +++ b/tests/cxx/unit/integrals/testing/testing.hpp @@ -0,0 +1,82 @@ +/* Copyright 2024 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 +#include +#include +#include +#include + +#include "ao_bases.hpp" +#include "molecules.hpp" +#include "shell_quartets.hpp" + +namespace integrals::testing { + +inline auto initialize_integrals() { + auto mm = pluginplay::ModuleManager(); + integrals::load_modules(mm); + integrals::set_defaults(mm); + return mm; +} + +template +auto eigen_tensor_(const tensorwrapper::buffer::BufferBase& buffer, + std::array extents, std::index_sequence) { + using namespace tensorwrapper; + const auto pdata = buffer::get_raw_data(buffer); + using eigen_type = Eigen::Tensor; + return Eigen::TensorMap(pdata.data(), extents[Is]...); +} + +// Checking eigen outputs +template +auto eigen_tensor(const tensorwrapper::buffer::BufferBase& buffer) { + std::array extents; + auto shape = buffer.layout().shape().as_smooth(); + for(std::size_t i = 0; i < N; ++i) extents[i] = shape.extent(i); + return eigen_tensor_(buffer, extents, + std::make_index_sequence()); +} + +template +auto trace(const tensorwrapper::buffer::BufferBase& buffer) { + auto t = eigen_tensor(buffer); + Eigen::Tensor trace = t.trace(); + return trace.coeff(); +} + +template +auto norm(const tensorwrapper::buffer::BufferBase& buffer) { + auto t = eigen_tensor(buffer); + Eigen::Tensor norm = t.square().sum().sqrt(); + return norm.coeff(); +} + +inline auto h2_mos() { + using mos_type = simde::type::mos; + using tensor_type = typename mos_type::transform_type; + tensor_type c({{-0.565516, -1.07019}, {-0.565516, 1.07019}}); + return mos_type(simde::type::aos(h2_sto3g_basis_set()), std::move(c)); +} + +inline auto h2_density() { + using density_type = simde::type::decomposable_e_density; + typename density_type::value_type rho( + {{0.31980835, 0.31980835}, {0.31980835, 0.31980835}}); + return density_type(rho, h2_mos()); +} + +} // namespace integrals::testing \ No newline at end of file diff --git a/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp b/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp index d86f81c8..07657803 100644 --- a/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp +++ b/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp @@ -13,10 +13,10 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#include "../libint/test_error.hpp" +#include "../testing/testing.hpp" #include -using namespace integrals::libint::test; +using namespace integrals::testing; using namespace integrals::utils; using eri_pt = simde::ERI4; diff --git a/tests/cxx/unit/integrals/utils/test_decontract_basis_set.cpp b/tests/cxx/unit/integrals/utils/test_decontract_basis_set.cpp index 774617d6..c936ad4f 100644 --- a/tests/cxx/unit/integrals/utils/test_decontract_basis_set.cpp +++ b/tests/cxx/unit/integrals/utils/test_decontract_basis_set.cpp @@ -14,18 +14,19 @@ * limitations under the License. */ -#include "../testing.hpp" +#include "../testing/testing.hpp" #include +using namespace integrals::testing; + TEST_CASE("DecontractBasisSet") { using test_pt = integrals::property_types::DecontractBasisSet; - pluginplay::ModuleManager mm; - integrals::load_modules(mm); + auto mm = initialize_integrals(); REQUIRE(mm.count("Decontract Basis Set")); // Get basis set - auto aobs = test::water_sto3g_basis_set(); + auto aobs = water_sto3g_basis_set(); // Get module auto& mod = mm.at("Decontract Basis Set"); @@ -34,6 +35,6 @@ TEST_CASE("DecontractBasisSet") { auto decontracted_aobs = mod.run_as(aobs); // Check output - auto corr = test::water_decontracted_sto3g_basis_set(); + auto corr = water_decontracted_sto3g_basis_set(); REQUIRE(decontracted_aobs == corr); } From 2decb239fff4b09cfe5ec1f6f144bf746de165ee Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Mon, 23 Feb 2026 11:02:02 -0600 Subject: [PATCH 06/14] cs screening seems to work --- .../libint/black_box_primitive_estimator.cpp | 3 +- .../cauchy_schwarz_primitive_estimator.cpp | 139 +++++++++++------- src/integrals/utils/rank2_shell_norm.hpp | 2 +- .../libint/black_box_primitive_estimator.cpp | 5 + .../cauchy_schwarz_primitive_estimator.cpp | 137 +++++++++++++---- 5 files changed, 203 insertions(+), 83 deletions(-) diff --git a/src/integrals/libint/black_box_primitive_estimator.cpp b/src/integrals/libint/black_box_primitive_estimator.cpp index e2b85567..fc5172f5 100644 --- a/src/integrals/libint/black_box_primitive_estimator.cpp +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -97,7 +97,8 @@ MODULE_RUN(BlackBoxPrimitiveEstimator) { for(shells[0] = 0; shells[0] < n_shells[0]; ++shells[0]) { const auto& bra_shell = bra_libint.at(shells[0]); assert(bra_shell.contr.size() == 1); // No general contraction support - const auto& bra_coeff = bra_shell.contr[0].coeff; + const auto& bra_coeff = bra_shell.contr[0].coeff; + const auto& bra_alpha = bra_shell.alpha; const auto n_prims_bra_shell = bra_coeff.size(); diff --git a/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp index a312d947..c6665ca0 100644 --- a/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp +++ b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -15,6 +15,7 @@ */ #include "../utils/rank2_shell_norm.hpp" +#include "detail_/make_libint_basis_set.hpp" #include "libint.hpp" #include #include @@ -41,71 +42,105 @@ MODULE_CTOR(CauchySchwarzPrimitiveEstimator) { 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 = bra_prims.n_primitives(); - const auto n_ket = ket_prims.n_primitives(); + 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 == bra_basis.n_primitives()); - assert(n_ket == ket_basis.n_primitives()); + assert(n_bra_prims == bra_basis.n_primitives()); + assert(n_ket_prims == ket_basis.n_primitives()); - // TODO: Only get the hyper diagonal we actually need + // TODO: We only need the hyper diagonal, so this is very wasteful simde::type::aos_squared bra(bra_prims, ket_prims); simde::type::v_ee_type v_ee{}; chemist::braket::BraKet mnls(bra, v_ee, bra); const auto& prim4 = submods.at("ERI4").run_as(mnls); + // TODO: Make our basis set normalize itself. + 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 buffer(n_bra * n_ket, 0.0); - - using iter_type = std::decay_t; // Type of the loop index - - // We loop over the decontracted basis, but make sure to grab the - // contraction coefficients from the real basis set. By looping over the - // decontracted basis we can make this look like the "usual" shell pair - // loop - // N.b. there's one shell per primitive in the real basis - // N.b. since the basis is decontracted the number of Cartesian/spherical - // AOs in a shell is the same as the number of Cartesian/spherical - // primitives in that shell - - std::array offsets{0, 0, 0, 0}; - std::array naos{0, 0, 0, 0}; - for(iter_type shell_i = 0; shell_i < bra_prims.n_shells(); ++shell_i) { - const auto bra_prim_i = bra_basis.primitive(shell_i); - const auto bra_coeff = std::fabs(bra_prim_i.coefficient()); - naos[0] = bra_prims.shell(shell_i).size(); - naos[2] = naos[0]; - - const auto bra_shell_offset = shell_i * n_ket; - offsets[1] = 0; - offsets[3] = 0; - for(iter_type shell_j = 0; shell_j < ket_prims.n_shells(); ++shell_j) { - const auto ket_prim_i = ket_basis.primitive(shell_j); - const auto ket_coeff = std::fabs(ket_prim_i.coefficient()); - naos[1] = ket_prims.shell(shell_j).size(); - naos[3] = naos[1]; - - auto C_ab = bra_coeff * ket_coeff; - auto shell_norm = utils::rank2_shell_norm(eris, offsets, naos); - buffer[bra_shell_offset + shell_j] = C_ab * std::sqrt(shell_norm); - - offsets[1] += naos[1]; - offsets[3] += naos[3]; - } - offsets[0] += naos[0]; - offsets[2] += naos[2]; - } - - tensorwrapper::shape::Smooth shape({n_bra, n_ket}); - tensorwrapper::buffer::Contiguous tw_buffer(std::move(buffer), shape); - simde::type::tensor rv(shape, std::move(tw_buffer)); + 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 + + 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]; + + 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]; + + } // loop over bra primitives + + prim_offsets[0] += n_prims_bra_shell; + + } // loop over bra shells + + simde::type::tensor rv(shape, std::move(buffer)); auto result = results(); return pt::wrap_results(result, rv); diff --git a/src/integrals/utils/rank2_shell_norm.hpp b/src/integrals/utils/rank2_shell_norm.hpp index ab522569..e97fa252 100644 --- a/src/integrals/utils/rank2_shell_norm.hpp +++ b/src/integrals/utils/rank2_shell_norm.hpp @@ -10,7 +10,7 @@ template auto rank2_shell_norm(TensorType&& buffer, std::array offsets, std::array naos) { using float_type = double; - // Take the shell norm over the primitive integrals in this shell pair. + float_type shell_norm = 0.0; std::array ao{0, 0, 0, 0}; std::vector abs_ao(4, 0); // Absolute indices diff --git a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp index 52134ed5..c6a52b74 100644 --- a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp +++ b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp @@ -23,6 +23,11 @@ using pt = integrals::property_types::PrimitivePairEstimator; namespace { +/* Correct values come from a MWE I made that recreates the integral screening + * using libint directly. Since the MWE exactly matches libint, I strongly + * believe these values are also correct. + */ + template auto corr_k01() { std::vector buffer{3.693e-38, 4.76823e-13, 8.9425e-06, diff --git a/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp index 252b1e4c..9f297a5b 100644 --- a/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp +++ b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -1,29 +1,108 @@ -// /* -// * 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 "../testing/testing.hpp" -// #include - -// using pt = integrals::property_types::PrimitivePairEstimator; -// using namespace integrals::testing; -// TEST_CASE("CauchySchwarzPrimitiveEstimator") { -// auto mm = initialize_integrals(); -// auto& mod = mm.at("CauchySchwarz Estimator"); -// auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); - -// auto Q_ab = mod.run_as(bra0, bra1); -// auto Q_cd = mod.run_as(ket0, ket1); -// // std::cout << Q_ab << std::endl; -// } \ No newline at end of file +/* + * 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 "../testing/testing.hpp" +#include + +using pt = integrals::property_types::PrimitivePairEstimator; +using namespace integrals::testing; +using tensorwrapper::operations::approximately_equal; + +namespace { + +/* Corr values come from a companion Python script that uses the dumped ERIs. + * The script uses Numpy and is much cleaner than the C++ code, which decreases + * the likelihood of errors in the corr values. That said, it is possible that + * there is an error in how I approached this. If such an error occurred it + * would likely show up in both places and NOT be caught by these tests. + */ + +// Assumes all s-type functions +template +auto Q_ab_corr() { + std::vector data{0.00000000e+00, 0.00000000e+00, 3.96602449e-10, + 0.00000000e+00, 3.75905705e-15, 2.01602647e-08, + 3.96602449e-10, 2.01602647e-08, 8.48436812e-07}; + tensorwrapper::shape::Smooth shape({3, 3}); + tensorwrapper::buffer::Contiguous cont(std::move(data), shape); + return simde::type::tensor(std::move(cont), shape); +} + +// Assumes p for basis 0 and s for basis 1 functions +template +auto Q_ab_psps_corr() { + std::vector data{0.00000000e+00, 0.00000000e+00, 2.88455917e-09, + 0.00000000e+00, 1.93317006e-13, 1.96314576e-07, + 1.03768892e-08, 3.60782670e-07, 6.23308753e-06}; + tensorwrapper::shape::Smooth shape({3, 3}); + tensorwrapper::buffer::Contiguous cont(std::move(data), shape); + return simde::type::tensor(std::move(cont), shape); +} + +// Assumes all p-type functions +template +auto Q_ab_pppp_corr() { + std::vector data{0.00000000e+00, 0.00000000e+00, 6.85302613e-08, + 0.00000000e+00, 9.31619988e-12, 3.04372388e-06, + 6.85302613e-08, 3.04372388e-06, 3.65812603e-05}; + tensorwrapper::shape::Smooth shape({3, 3}); + tensorwrapper::buffer::Contiguous cont(std::move(data), shape); + return simde::type::tensor(std::move(cont), shape); +} + +template +auto Q_cd_corr() { + std::vector data{0.00000000, 2.08398137e-08, 3.10753955e-05, + 2.08398137e-08, 1.15521140e-05, 2.21832728e-04, + 3.10753955e-05, 2.21832728e-04, 3.13532082e-04}; + tensorwrapper::shape::Smooth shape({3, 3}); + tensorwrapper::buffer::Contiguous cont(std::move(data), shape); + return simde::type::tensor(std::move(cont), shape); +} + +} // namespace + +TEST_CASE("CauchySchwarzPrimitiveEstimator") { + auto mm = initialize_integrals(); + auto& mod = mm.at("CauchySchwarz Estimator"); + auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + + SECTION("Q(bra0, bra1)") { + auto Q_ab = mod.run_as(bra0, bra1); + auto corr = Q_ab_corr(); + REQUIRE(approximately_equal(Q_ab, corr, 1E-8)); + } + + SECTION("Q(ket0, ket1)") { + auto Q_cd = mod.run_as(ket0, ket1); + auto corr = Q_cd_corr(); + REQUIRE(approximately_equal(Q_cd, corr, 1E-8)); + } + + SECTION("Q_ab (ps|ps)") { + bra0.shell(0).l() = 1; + auto Q_ab = mod.run_as(bra0, bra1); + auto corr = Q_ab_psps_corr(); + REQUIRE(approximately_equal(Q_ab, corr, 1E-8)); + } + + SECTION("Q_ab (pp|pp)") { + bra0.shell(0).l() = 1; + bra1.shell(0).l() = 1; + auto Q_ab = mod.run_as(bra0, bra1); + auto corr = Q_ab_pppp_corr(); + REQUIRE(approximately_equal(Q_ab, corr, 1E-8)); + } +} \ No newline at end of file From 374130f6a7b6066873eb48fd25531686cd3b632d Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Mon, 23 Feb 2026 11:17:19 -0600 Subject: [PATCH 07/14] backup --- .../libint/primitive_error_model.cpp | 2 +- .../libint/primitive_error_model.cpp | 73 ++++++------------- 2 files changed, 24 insertions(+), 51 deletions(-) diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp index 873c11e1..e687ee68 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -67,7 +67,7 @@ auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, } // End loop over prim[2] } // End loop over prim[1] } // End loop over prim[0] - + std::cout << "Shell block error: " << error << std::endl; return error; } diff --git a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp index 4df662db..85edcdc6 100644 --- a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -1,50 +1,23 @@ -// #include "../testing/testing.hpp" -// #include - -// using eri4_pt = simde::ERI4; -// using pt = integrals::property_types::Uncertainty; - -// using namespace integrals::testing; - -// TEST_CASE("PrimitiveErrorModel") { -// auto mm = initialize_integrals(); - -// auto& mod = mm.at("Primitive Error Model"); - -// simde::type::v_ee_type v_ee{}; -// using float_type = double; -// SECTION("H2") { -// auto h2_aos = test::h2_sto3g_basis_set(); -// simde::type::aos_squared h2_aos2(h2_aos, h2_aos); -// chemist::braket::BraKet mnls(h2_aos2, v_ee, h2_aos2); -// auto error = mod.run_as(mnls, 1.0e-10); - -// tensorwrapper::shape::Smooth shape({2, 2, 2, 2}); -// std::vector data(shape.size()); -// tensorwrapper::buffer::Contiguous buffer(std::move(data), shape); -// REQUIRE(buffer.approximately_equal(error.buffer(), 1.0e-12)); -// } - -// SECTION("H2 2") { -// float_type tol = 1.0e-10; -// auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); -// auto value = manual_contract_shell(bra0, bra1, ket0, ket1, mm); -// std::cout << "Screened value: " << value << std::endl; - -// simde::type::aos_squared bra(bra0, bra1); -// simde::type::aos_squared ket(ket0, ket1); -// chemist::braket::BraKet mnls(bra, v_ee, ket); - -// auto error = mod.run_as(mnls, tol); -// auto corr = mm.at("Analytic Error").run_as(mnls, tol); - -// auto eri4_mod = mm.at("ERI4").unlocked_copy(); -// eri4_mod.change_input("Threshold", float_type(1E-16)); -// auto i_value = eri4_mod.run_as(mnls); -// std::cout << "Corr Screened Value: " << i_value << std::endl; - -// // Our value: 0.0002263495591894 -// // Libint's value: 0.0002263440759384 -// // No screening: 0.0002263495626484 -// } -// } \ No newline at end of file +#include "../testing/testing.hpp" +#include + +using eri4_pt = simde::ERI4; +using pt = integrals::property_types::Uncertainty; + +using namespace integrals::testing; + +TEST_CASE("PrimitiveErrorModel") { + auto mm = initialize_integrals(); + + auto& mod = mm.at("Primitive Error Model"); + auto& anal_error = mm.at("Analytic Error"); + auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + simde::type::aos_squared bra(bra0, bra1); + simde::type::v_ee_type v_ee{}; + simde::type::aos_squared ket(ket0, ket1); + chemist::braket::BraKet mnls(bra, v_ee, ket); + double tol = 1E-10; + auto error = mod.run_as(mnls, tol); + std::cout << error << std::endl; + std::cout << anal_error.run_as(mnls, tol) << std::endl; +} \ No newline at end of file From 8814ac14096380cb4a3111df7d58658c88617411 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Mon, 23 Feb 2026 13:19:09 -0600 Subject: [PATCH 08/14] module for screening pairs --- include/integrals/property_types.hpp | 18 ++++++ src/integrals/ao_integrals/ao_integrals.hpp | 2 +- src/integrals/libint/libint.cpp | 6 +- .../libint/primitive_error_model.cpp | 1 - .../utils/screen_primitive_pairs.cpp | 58 +++++++++++++++++++ src/integrals/utils/utils.hpp | 7 ++- .../integrals/ao_integrals/test_uq_driver.cpp | 16 +++-- .../libint/primitive_error_model.cpp | 6 +- .../utils/screen_primitive_pairs.cpp | 37 ++++++++++++ 9 files changed, 138 insertions(+), 13 deletions(-) create mode 100644 src/integrals/utils/screen_primitive_pairs.cpp create mode 100644 tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp diff --git a/include/integrals/property_types.hpp b/include/integrals/property_types.hpp index 2798c44f..3998de81 100644 --- a/include/integrals/property_types.hpp +++ b/include/integrals/property_types.hpp @@ -54,6 +54,24 @@ PROPERTY_TYPE_RESULTS(PrimitivePairEstimator) { return rv; } +DECLARE_PROPERTY_TYPE(PairScreener); +PROPERTY_TYPE_INPUTS(PairScreener) { + using ao_basis = const simde::type::ao_basis_set&; + auto rv = pluginplay::declare_input() + .add_field("Bra Basis Set") + .add_field("Ket Basis Set") + .add_field("Tolerance"); + return rv; +} + +PROPERTY_TYPE_RESULTS(PairScreener) { + using index_vector = std::vector; + using pair_vector = std::vector; + auto rv = pluginplay::declare_result().add_field( + "Primitive Pairs Passing Screening"); + return rv; +} + template DECLARE_TEMPLATED_PROPERTY_TYPE(Uncertainty, BasePT); diff --git a/src/integrals/ao_integrals/ao_integrals.hpp b/src/integrals/ao_integrals/ao_integrals.hpp index cee09ef4..6fe9fb10 100644 --- a/src/integrals/ao_integrals/ao_integrals.hpp +++ b/src/integrals/ao_integrals/ao_integrals.hpp @@ -40,7 +40,7 @@ inline void set_defaults(pluginplay::ModuleManager& mm) { mm.change_submod("Density Fitting Integral", "Coulomb Metric", "Coulomb Metric"); mm.change_submod("UQ Driver", "ERIs", "ERI4"); - mm.change_submod("UQ Driver", "ERI Error", "Analytic Error"); + mm.change_submod("UQ Driver", "ERI Error", "Primitive Error Model"); } inline void load_modules(pluginplay::ModuleManager& mm) { diff --git a/src/integrals/libint/libint.cpp b/src/integrals/libint/libint.cpp index b63618e8..952ecdc2 100644 --- a/src/integrals/libint/libint.cpp +++ b/src/integrals/libint/libint.cpp @@ -199,8 +199,10 @@ void set_defaults(pluginplay::ModuleManager& mm) { "CauchySchwarz Estimator"); mm.change_submod("CauchySchwarz Estimator", "Decontract Basis Set", "Decontract Basis Set"); - mm.change_submod("CauchySchwarz Estimator", "ERI4", "ERI4"); - mm.change_submod("Analytic Error", "ERI4s", "ERI4"); + mm.copy_module("ERI4", "Benchmark ERI4"); + mm.change_input("Benchmark ERI4", "Threshold", 1.0E-16); + mm.change_submod("CauchySchwarz Estimator", "ERI4", "Benchmark ERI4"); + mm.change_submod("Analytic Error", "ERI4s", "Benchmark ERI4"); } #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 e687ee68..b7be89ac 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -67,7 +67,6 @@ auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, } // End loop over prim[2] } // End loop over prim[1] } // End loop over prim[0] - std::cout << "Shell block error: " << error << std::endl; return error; } diff --git a/src/integrals/utils/screen_primitive_pairs.cpp b/src/integrals/utils/screen_primitive_pairs.cpp new file mode 100644 index 00000000..a4cc06c8 --- /dev/null +++ b/src/integrals/utils/screen_primitive_pairs.cpp @@ -0,0 +1,58 @@ +#include "utils.hpp" +#include +namespace integrals::utils { + +namespace { + +const auto desc = R"( +Screen Primitive Pairs +====================== + +This module returns a set of primitive pairs whose estimated contribution to +the overall ERI is greater than a user specified threshold. +)"; + +} + +using pair_estimate_pt = integrals::property_types::PrimitivePairEstimator; +using pt = integrals::property_types::PairScreener; + +MODULE_CTOR(ScreenPrimitivePairs) { + description(desc); + satisfies_property_type(); + + add_submodule("Primitive Pair Estimator") + .set_description("The module used to estimate the contributions of " + "primitive pairs to the overall integral values"); +} + +MODULE_RUN(ScreenPrimitivePairs) { + const auto&& [bra, ket, tol] = pt::unwrap_inputs(inputs); + + auto& screener = submods.at("Primitive Pair Estimator"); + const auto& Kij = screener.run_as(bra, ket); + + const auto n_bra_prims = bra.n_primitives(); + const auto n_ket_prims = ket.n_primitives(); + + using float_type = double; + using tensorwrapper::buffer::make_contiguous; + const auto& buffer = make_contiguous(Kij.buffer()); + + using index_vector = std::vector; + index_vector prim{0, 0}; + std::vector rv; + using wtf::fp::float_cast; + + for(prim[0] = 0; prim[0] < n_bra_prims; ++prim[0]) { + for(prim[1] = 0; prim[1] < n_ket_prims; ++prim[1]) { + const auto val = float_cast(buffer.get_elem(prim)); + if(val > tol) rv.push_back(prim); + } + } + + auto result = results(); + return pt::wrap_results(result, rv); +} + +} // namespace integrals::utils \ No newline at end of file diff --git a/src/integrals/utils/utils.hpp b/src/integrals/utils/utils.hpp index 5109d29e..72e3cc36 100644 --- a/src/integrals/utils/utils.hpp +++ b/src/integrals/utils/utils.hpp @@ -30,11 +30,16 @@ namespace integrals::utils { DECLARE_MODULE(DecontractBasisSet); +DECLARE_MODULE(ScreenPrimitivePairs); -inline void set_defaults(pluginplay::ModuleManager& mm) {} +inline void set_defaults(pluginplay::ModuleManager& mm) { + mm.change_submod("Screen Primitive Pairs", "Primitive Pair Estimator", + "Black Box Primitive Pair Estimator"); +} inline void load_modules(pluginplay::ModuleManager& mm) { mm.add_module("Decontract Basis Set"); + mm.add_module("Screen Primitive Pairs"); set_defaults(mm); } diff --git a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp index 398d6cec..51da1568 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp @@ -57,7 +57,8 @@ TEST_CASE("UQ Driver") { integrals::load_modules(mm); integrals::set_defaults(mm); REQUIRE(mm.count("UQ Driver")); - + mm.copy_module("UQ Driver", "UQ w/analytic Error"); + mm.change_submod("UQ w/analytic Error", "ERI Error", "Analytic Error"); mm.change_input("ERI4", "Threshold", 1.0e-6); // Get basis set @@ -75,9 +76,14 @@ TEST_CASE("UQ Driver") { chemist::braket::BraKet braket(aos_squared, op, aos_squared); // Call module - auto T = mm.at("UQ Driver").run_as(braket); - auto T_corr = corr_answer(T); - using tensorwrapper::operations::approximately_equal; - REQUIRE(approximately_equal(T_corr, T, 1E-6)); + // auto T = mm.at("UQ Driver").run_as(braket); + // auto T_corr = mm.at("UQ w/analytic Error").run_as(braket); + + // std::cout << T << std::endl; + // std::cout << T_corr << std::endl; + + // auto T_corr = corr_answer(T); + // using tensorwrapper::operations::approximately_equal; + // REQUIRE(approximately_equal(T_corr, T, 1E-6)); } } diff --git a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp index 85edcdc6..b4291d48 100644 --- a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -17,7 +17,7 @@ TEST_CASE("PrimitiveErrorModel") { simde::type::aos_squared ket(ket0, ket1); chemist::braket::BraKet mnls(bra, v_ee, ket); double tol = 1E-10; - auto error = mod.run_as(mnls, tol); - std::cout << error << std::endl; - std::cout << anal_error.run_as(mnls, tol) << std::endl; + // auto error = mod.run_as(mnls, tol); + // std::cout << error << std::endl; + // std::cout << anal_error.run_as(mnls, tol) << std::endl; } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp b/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp new file mode 100644 index 00000000..a12747c7 --- /dev/null +++ b/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp @@ -0,0 +1,37 @@ +/* + * 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 "../testing/testing.hpp" +#include + +using namespace integrals::testing; + +using pt = integrals::property_types::PairScreener; + +TEST_CASE("Screen Primitive Pairs") { + auto mm = initialize_integrals(); + + auto mod = mm.at("Screen Primitive Pairs"); + + auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + + SECTION("(Bra0, Bra1) pairs") { + auto pairs = mod.run_as(bra0, bra1, 1E-10); + for(const auto& pair : pairs) { + for(const auto& prim : pair) { std::cout << prim << " "; } + std::cout << std::endl; + } + } +} \ No newline at end of file From e6cebb0090f1e6cc758ffea9c378f81720cf4dda Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 24 Feb 2026 10:48:38 -0600 Subject: [PATCH 09/14] finally r2g --- CMakeLists.txt | 6 ++ examples/primitive_error_models/main.cpp | 74 +++++++++++++ .../libint/black_box_primitive_estimator.cpp | 11 +- src/integrals/libint/libint.cpp | 2 - .../libint/primitive_error_model.cpp | 29 ++--- .../integrals/ao_integrals/test_uq_driver.cpp | 6 +- .../unit/integrals/libint/analytic_error.cpp | 17 ++- .../libint/black_box_primitive_estimator.cpp | 73 ++++++++++--- .../libint/primitive_error_model.cpp | 100 +++++++++++++++--- .../unit/integrals/testing/shell_quartets.hpp | 53 ++++++++++ .../utils/screen_primitive_pairs.cpp | 40 ++++++- 11 files changed, 350 insertions(+), 61 deletions(-) create mode 100644 examples/primitive_error_models/main.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index eab806b6..95bc4b37 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -72,6 +72,12 @@ cmaize_add_library( DEPENDS "${project_depends}" ) +cmaize_add_executable( + primitive_error_models + SOURCE_DIR "examples/primitive_error_models" + DEPENDS "${PROJECT_NAME}" +) + include(nwx_pybind11) nwx_add_pybind11_module( ${PROJECT_NAME} diff --git a/examples/primitive_error_models/main.cpp b/examples/primitive_error_models/main.cpp new file mode 100644 index 00000000..bd2aa627 --- /dev/null +++ b/examples/primitive_error_models/main.cpp @@ -0,0 +1,74 @@ +#include +#include + +/* This example showcases how to: + * + * 1. Compute the analytic error in an ERI4 integral tensor owing to primitive + * pair screening. + */ + +namespace { + +// This makes a basis set for H2 (bond distance 1.40 a.u.) using STO-3G. +inline simde::type::ao_basis_set h2_sto3g_basis_set() { + using ao_basis_t = simde::type::ao_basis_set; + using atomic_basis_t = simde::type::atomic_basis_set; + using cg_t = simde::type::contracted_gaussian; + using point_t = simde::type::point; + using doubles_t = std::vector; + + point_t r0{0.0, 0.0, 0.0}; + point_t r1{0.0, 0.0, 1.40}; + + doubles_t cs{0.1543289673, 0.5353281423, 0.4446345422}; + doubles_t es{3.425250914, 0.6239137298, 0.1688554040}; + cg_t cg0(cs.begin(), cs.end(), es.begin(), es.end(), r0); + cg_t cg1(cs.begin(), cs.end(), es.begin(), es.end(), r1); + atomic_basis_t h0("sto-3g", 1, r0); + atomic_basis_t h1("sto-3g", 1, r1); + h0.add_shell(chemist::ShellType::cartesian, 0, cg0); + h1.add_shell(chemist::ShellType::cartesian, 0, cg1); + + ao_basis_t bs; + bs.add_center(h0); + bs.add_center(h1); + return bs; +} + +} // namespace + +// Property types for the ERI4 and the error in the ERI4 +using eri4_pt = simde::ERI4; +using eri4_error_pt = integrals::property_types::Uncertainty; + +int main(int argc, char* argv[]) { + // Makes sure the environment doesn't go out of scope before the end. + auto rt = std::make_unique(); + + // Initializes a ModuleManager object with the integrals plugin + pluginplay::ModuleManager mm(std::move(rt), nullptr); + integrals::load_modules(mm); + integrals::set_defaults(mm); + + // Modules for computing analytic error and estimating error + auto& analytic_error_mod = mm.at("Analytic Error"); + auto& error_model = mm.at("Primitive Error Model"); + + // Makes: basis set, direct product of the basis set, and 1/r12 operator + simde::type::aos aos(h2_sto3g_basis_set()); + simde::type::aos_squared aos2(aos, aos); + simde::type::v_ee_type op{}; + + // Make BraKet + chemist::braket::BraKet mnls(aos2, op, aos2); + + // Compute the error by screening with tolerance "tol" + double tol = 1E-10; + auto error = analytic_error_mod.run_as(mnls, tol); + auto approx_error = error_model.run_as(mnls, tol); + + std::cout << "Analytic error: " << error << std::endl; + std::cout << "Estimated error: " << approx_error << std::endl; + + return 0; +} \ No newline at end of file diff --git a/src/integrals/libint/black_box_primitive_estimator.cpp b/src/integrals/libint/black_box_primitive_estimator.cpp index fc5172f5..2c76c25d 100644 --- a/src/integrals/libint/black_box_primitive_estimator.cpp +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -59,6 +59,14 @@ auto compute_k(T zeta_i, T zeta_j, T coeff_i, T coeff_j, T dr2) { return coeff_i * coeff_j * std::exp(ratio * dr2); } +template +auto compute_ln_k(T zeta_i, T zeta_j, T coeff_i, T coeff_j, T dr2) { + const auto num = -zeta_i * zeta_j; + const auto denom = zeta_i + zeta_j; + const auto ratio = num / denom; + return std::log(coeff_i) + std::log(coeff_j) + ratio * dr2; +} + } // namespace using pt = integrals::property_types::PrimitivePairEstimator; @@ -97,8 +105,7 @@ MODULE_RUN(BlackBoxPrimitiveEstimator) { for(shells[0] = 0; shells[0] < n_shells[0]; ++shells[0]) { const auto& bra_shell = bra_libint.at(shells[0]); assert(bra_shell.contr.size() == 1); // No general contraction support - const auto& bra_coeff = bra_shell.contr[0].coeff; - + const auto& bra_coeff = bra_shell.contr[0].coeff; const auto& bra_alpha = bra_shell.alpha; const auto n_prims_bra_shell = bra_coeff.size(); diff --git a/src/integrals/libint/libint.cpp b/src/integrals/libint/libint.cpp index 952ecdc2..57fadc32 100644 --- a/src/integrals/libint/libint.cpp +++ b/src/integrals/libint/libint.cpp @@ -195,8 +195,6 @@ void set_defaults(pluginplay::ModuleManager& mm) { mm.change_submod("Primitive Error Model", "Black Box Primitive Pair Estimator", "Black Box Primitive Pair Estimator"); - mm.change_submod("Primitive Error Model", "Primitive Pair Estimator", - "CauchySchwarz Estimator"); mm.change_submod("CauchySchwarz Estimator", "Decontract Basis Set", "Decontract Basis Set"); mm.copy_module("ERI4", "Benchmark ERI4"); diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp index b7be89ac..6cdce089 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -10,8 +10,7 @@ const auto desc = "Compute uncertainty estimates for ERI4 integrals"; // Sums contributions Q_abQ_cd when Q_ab or Q_cd is below tol template auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, - TensorType&& Q_ab, TensorType&& Q_cd, TensorType&& K_ab, - TensorType&& K_cd, double tol) { + TensorType&& K_ab, TensorType&& K_cd, double tol) { // Get the number of primitives in each shell const auto n_prim0 = shells[0].n_primitives(); const auto n_prim1 = shells[1].n_primitives(); @@ -40,8 +39,6 @@ auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, auto K_01 = float_cast(K_ab.get_elem(i01)); auto K_01_mag = std::fabs(K_01); bool K_01_neglected = K_01_mag < tol; - auto Q_01 = float_cast(Q_ab.get_elem(i01)); - auto Q_01_mag = std::fabs(Q_01); for(prim[2] = 0; prim[2] < n_prim2; ++prim[2]) { const auto p2 = prim_offsets[2] + prim[2]; @@ -54,14 +51,12 @@ auto shell_block_error(ShellsType&& shells, OffsetTypes&& prim_offsets, auto K_23 = float_cast(K_cd.get_elem(i23)); auto K_23_mag = std::fabs(K_23); bool K_23_neglected = K_23_mag < tol; - auto Q_23 = float_cast(Q_cd.get_elem(i23)); - auto Q_23_mag = std::fabs(Q_23); auto prod_neglected = K_01_mag * K_23_mag < tol; // If either was neglected we pick up an error Q_01 * Q_23 if(K_01_neglected || K_23_neglected || prod_neglected) { - error += Q_01_mag * Q_23_mag; + error += tol; } } // End loop over prim[3] } // End loop over prim[2] @@ -118,9 +113,6 @@ MODULE_CTOR(PrimitiveErrorModel) { add_submodule("Black Box Primitive Pair Estimator") .set_description("The module used to estimate the contributions of " "primitive pairs to the overall integral values"); - add_submodule("Primitive Pair Estimator") - .set_description("The module used to estimate the contributions of " - "primitive pairs to the overall integral values"); } MODULE_RUN(PrimitiveErrorModel) { @@ -131,21 +123,14 @@ MODULE_RUN(PrimitiveErrorModel) { const auto ket0 = braket.ket().first.ao_basis_set(); const auto ket1 = braket.ket().second.ao_basis_set(); - // Get the Q_ab and Q_cd tensors (used to estimate error) - auto& estimator = submods.at("Primitive Pair Estimator"); - const auto& Q_ab_tw = estimator.run_as(bra0, bra1); - const auto& Q_cd_tw = estimator.run_as(ket0, ket1); - // Get the K_ab and K_cd tensors (used to determine which primitives are // neglected) - auto& estimator2 = submods.at("Black Box Primitive Pair Estimator"); - const auto& K_ab_tw = estimator2.run_as(bra0, bra1); - const auto& K_cd_tw = estimator2.run_as(ket0, ket1); + auto& estimator = submods.at("Black Box Primitive Pair Estimator"); + const auto& K_ab_tw = estimator.run_as(bra0, bra1); + const auto& K_cd_tw = estimator.run_as(ket0, ket1); // XXX: Workaround for needing the contiguous buffers to access elements using tensorwrapper::buffer::make_contiguous; - const auto& Q_ab = make_contiguous(Q_ab_tw.buffer()); - const auto& Q_cd = make_contiguous(Q_cd_tw.buffer()); const auto& K_ab = make_contiguous(K_ab_tw.buffer()); const auto& K_cd = make_contiguous(K_cd_tw.buffer()); @@ -191,8 +176,8 @@ MODULE_RUN(PrimitiveErrorModel) { for(shell_i[3] = 0; shell_i[3] < n_shells[3]; ++shell_i[3]) { shells[3] = ket1.shell(shell_i[3]); - auto shell_error = shell_block_error( - shells, prim_offset, Q_ab, Q_cd, K_ab, K_cd, tol); + auto shell_error = + shell_block_error(shells, prim_offset, K_ab, K_cd, tol); fill_ao_block(shells, ao_offsets, buffer, shell_error); prim_offset[3] += shells[3].n_primitives(); diff --git a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp index 51da1568..e85ab324 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_uq_driver.cpp @@ -75,9 +75,9 @@ TEST_CASE("UQ Driver") { // Make BraKet Input chemist::braket::BraKet braket(aos_squared, op, aos_squared); - // Call module - // auto T = mm.at("UQ Driver").run_as(braket); - // auto T_corr = mm.at("UQ w/analytic Error").run_as(braket); + // Call modules + auto T = mm.at("UQ Driver").run_as(braket); + auto T_corr = mm.at("UQ w/analytic Error").run_as(braket); // std::cout << T << std::endl; // std::cout << T_corr << std::endl; diff --git a/tests/cxx/unit/integrals/libint/analytic_error.cpp b/tests/cxx/unit/integrals/libint/analytic_error.cpp index 1caaa132..071d9889 100644 --- a/tests/cxx/unit/integrals/libint/analytic_error.cpp +++ b/tests/cxx/unit/integrals/libint/analytic_error.cpp @@ -9,8 +9,7 @@ using tensorwrapper::operations::approximately_equal; namespace { template -auto corr_error() { - FloatType diff = -0.0000000054867100; +auto corr_error(FloatType diff) { tensorwrapper::shape::Smooth shape({1, 1, 1, 1}); std::vector buffer(shape.size(), diff); tensorwrapper::buffer::Contiguous cont(std::move(buffer), shape); @@ -34,7 +33,19 @@ TEST_CASE("AnalyticError") { chemist::braket::BraKet mnls(bra, v_ee, ket); float_type tol = 1.0e-10; auto error = mod.run_as(mnls, tol); - auto corr = corr_error(); + auto corr = corr_error(-0.0000000054867100); + REQUIRE(approximately_equal(error, corr, 1.0e-12)); + } + + SECTION("H2O/STO-3G (00|34)") { + auto [bra0, bra1, ket0, ket1] = get_h2o_0034_bases(); + + simde::type::aos_squared bra(bra0, bra1); + simde::type::aos_squared ket(ket0, ket1); + chemist::braket::BraKet mnls(bra, v_ee, ket); + float_type tol = 1.0e-10; + auto error = mod.run_as(mnls, tol); + auto corr = corr_error(-0.0000000000315610); REQUIRE(approximately_equal(error, corr, 1.0e-12)); } } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp index c6a52b74..3eb07127 100644 --- a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp +++ b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp @@ -23,7 +23,8 @@ using pt = integrals::property_types::PrimitivePairEstimator; namespace { -/* Correct values come from a MWE I made that recreates the integral screening +/* Correct values for H2 come from a MWE I made that recreates the integral + * screening * using libint directly. Since the MWE exactly matches libint, I strongly * believe these values are also correct. */ @@ -49,24 +50,70 @@ auto corr_k23() { return simde::type::tensor(std::move(cont), shape); } +/* Correct values for water comes from a companion python script that I am less + * confident in than the H2 values. The Python script is much simpler than the + * C++ implementation, and I didn't see any obvious errors, so it's quite + * plausible that it's correct, but if I made a logic error in the C++ code I + * likely made the same logic error int the Python code. + */ +template +auto corr_k01_water() { + std::vector buffer{ + 18.0790216813702180, 17.4852396535751247, 5.4493863186785507, + 17.4852396535751247, 16.9109596266485731, 5.2704082901341138, + 5.4493863186785507, 5.2704082901341138, 1.6425563160202101}; + tensorwrapper::shape::Smooth shape({3, 3}); + tensorwrapper::buffer::Contiguous cont(std::move(buffer), shape); + return simde::type::tensor(std::move(cont), shape); +} + +template +auto corr_k23_water() { + std::vector buffer{ + 0.0000000007980106, 0.0002571668526586, 0.0041100638983562, + 0.0002571668526586, 0.0025216231166686, 0.0053704017775550, + 0.0041100638983562, 0.0053704017775550, 0.0028156052053020}; + tensorwrapper::shape::Smooth shape({3, 3}); + tensorwrapper::buffer::Contiguous cont(std::move(buffer), shape); + return simde::type::tensor(std::move(cont), shape); +} + } // namespace TEST_CASE("BlackBoxPrimitiveEstimator") { - auto mm = initialize_integrals(); - auto& mod = mm.at("Black Box Primitive Pair Estimator"); - auto&& [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); - + auto mm = initialize_integrals(); + auto& mod = mm.at("Black Box Primitive Pair Estimator"); using float_type = double; - SECTION("K(bra0, bra1)") { - auto K01 = mod.run_as(bra0, bra1); - auto corr = corr_k01(); - REQUIRE(approximately_equal(K01, corr, 1E-8)); + SECTION("H2 Dimer/STO-3G (03|12)") { + auto&& [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + + SECTION("K(bra0, bra1)") { + auto K01 = mod.run_as(bra0, bra1); + auto corr = corr_k01(); + REQUIRE(approximately_equal(K01, corr, 1E-8)); + } + + SECTION("K(ket0, ket1)") { + auto K23 = mod.run_as(ket0, ket1); + auto corr = corr_k23(); + REQUIRE(approximately_equal(K23, corr, 1E-8)); + } } - SECTION("K(ket0, ket1)") { - auto K23 = mod.run_as(ket0, ket1); - auto corr = corr_k23(); - REQUIRE(approximately_equal(K23, corr, 1E-8)); + SECTION("H2O/STO-3G (00|34)") { + auto&& [bra0, bra1, ket0, ket1] = get_h2o_0034_bases(); + + SECTION("K(bra0, bra1)") { + auto K01 = mod.run_as(bra0, bra1); + auto corr = corr_k01_water(); + REQUIRE(approximately_equal(K01, corr, 1E-8)); + } + + SECTION("K(ket0, ket1)") { + auto K23 = mod.run_as(ket0, ket1); + auto corr = corr_k23_water(); + REQUIRE(approximately_equal(K23, corr, 1E-8)); + } } } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp index b4291d48..efd91718 100644 --- a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -6,18 +6,94 @@ using pt = integrals::property_types::Uncertainty; using namespace integrals::testing; -TEST_CASE("PrimitiveErrorModel") { - auto mm = initialize_integrals(); +namespace { + +template +bool compare_magnitudes(T0&& lhs, T1&& corr, double tol) { + using float_type = double; + using tensorwrapper::buffer::get_raw_data; + const auto lhs_buffer = get_raw_data(lhs.buffer()); + const auto corr_buffer = get_raw_data(corr.buffer()); + assert(lhs_buffer.size() == corr_buffer.size()); + + for(std::size_t i = 0; i < lhs_buffer.size(); ++i) { + auto lhs_val = lhs_buffer[i]; + auto corr_val = corr_buffer[i]; + if(lhs_val == 0.0 && corr_val == 0.0) { + continue; + } else if(lhs_val == 0.0 || corr_val == 0.0) { + auto abs_diff = std::fabs(lhs_val - corr_val); + if(abs_diff > 10 * tol) { + std::cout << "Absolute diff: " << abs_diff << std::endl; + return false; + } + continue; + } + auto lhs_mag = std::log10(std::fabs(lhs_val)); + auto corr_mag = std::log10(std::fabs(corr_val)); + if(lhs_mag != Catch::Approx(corr_mag).epsilon(0.2)) { + auto pct_error = (lhs_mag - corr_mag) / corr_mag; + std::cout << lhs_mag << " " << corr_mag << " " << pct_error + << std::endl; + return false; + } + } + return true; +} + +} // namespace - auto& mod = mm.at("Primitive Error Model"); - auto& anal_error = mm.at("Analytic Error"); - auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); - simde::type::aos_squared bra(bra0, bra1); +TEST_CASE("PrimitiveErrorModel") { + auto mm = initialize_integrals(); + auto& mod = mm.at("Primitive Error Model"); + auto& anal_error = mm.at("Analytic Error"); simde::type::v_ee_type v_ee{}; - simde::type::aos_squared ket(ket0, ket1); - chemist::braket::BraKet mnls(bra, v_ee, ket); - double tol = 1E-10; - // auto error = mod.run_as(mnls, tol); - // std::cout << error << std::endl; - // std::cout << anal_error.run_as(mnls, tol) << std::endl; + + SECTION("Single (ss|ss) quartet") { + auto [bra0, bra1, ket0, ket1] = get_h2_dimer_0312_bases(); + simde::type::aos_squared bra(bra0, bra1); + simde::type::aos_squared ket(ket0, ket1); + chemist::braket::BraKet mnls(bra, v_ee, ket); + + double tol = 1E-10; + auto error = mod.run_as(mnls, tol); + auto corr = anal_error.run_as(mnls, tol); + REQUIRE(compare_magnitudes(error, corr, tol)); + } + + SECTION("H2 molecule") { + 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); + + double tol = 1E-6; + auto error = mod.run_as(mnls, tol); + auto corr = anal_error.run_as(mnls, tol); + REQUIRE(compare_magnitudes(error, corr, tol)); + } + + SECTION("Water/STO-3G(00|34)") { + auto [bra0, bra1, ket0, ket1] = get_h2o_0034_bases(); + simde::type::aos_squared bra(bra0, bra1); + simde::type::aos_squared ket(ket0, ket1); + chemist::braket::BraKet mnls(bra, v_ee, ket); + + double tol = 1E-10; + auto error = mod.run_as(mnls, tol); + auto corr = anal_error.run_as(mnls, tol); + REQUIRE(compare_magnitudes(error, corr, tol)); + } + + SECTION("H2O molecule") { + auto aobs = water_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); + + double tol = 1E-10; + auto error = mod.run_as(mnls, tol); + auto corr = anal_error.run_as(mnls, tol); + REQUIRE(compare_magnitudes(error, corr, tol)); + } } \ No newline at end of file diff --git a/tests/cxx/unit/integrals/testing/shell_quartets.hpp b/tests/cxx/unit/integrals/testing/shell_quartets.hpp index 2670903e..48d99085 100644 --- a/tests/cxx/unit/integrals/testing/shell_quartets.hpp +++ b/tests/cxx/unit/integrals/testing/shell_quartets.hpp @@ -33,4 +33,57 @@ inline auto get_h2_dimer_0312_bases() { return std::make_tuple(bra0, bra1, ket0, ket1); } +/* After playing around I found that my code was getting hard zero for this + * quartet, but libint got a non-zero, but small value. + */ +inline auto get_h2o_0034_bases() { + simde::type::ao_basis_set bra0, ket3, ket4; + auto water_aos = water_sto3g_basis_set(); + + using abs_type = simde::type::ao_basis_set::value_type; + using cg_type = simde::type::contracted_gaussian; + using float_type = double; + using vector_type = std::vector; + + const auto shell0 = water_aos.shell(0); + const auto center0 = shell0.center().as_point(); + const auto cg0_view = shell0.contracted_gaussian(); + vector_type cs0{cg0_view[0].coefficient(), cg0_view[1].coefficient(), + cg0_view[2].coefficient()}; + vector_type es0{cg0_view[0].exponent(), cg0_view[1].exponent(), + cg0_view[2].exponent()}; + cg_type cg0(cs0.begin(), cs0.end(), es0.begin(), es0.end(), center0); + + abs_type abs0(center0); + abs0.add_shell(shell0.pure(), shell0.l(), cg0); + bra0.add_center(abs0); + + const auto shell3 = water_aos.shell(3); + const auto center3 = shell3.center().as_point(); + const auto cg3_view = shell3.contracted_gaussian(); + vector_type cs3{cg3_view[0].coefficient(), cg3_view[1].coefficient(), + cg3_view[2].coefficient()}; + vector_type es3{cg3_view[0].exponent(), cg3_view[1].exponent(), + cg3_view[2].exponent()}; + cg_type cg3(cs3.begin(), cs3.end(), es3.begin(), es3.end(), center3); + + abs_type abs3(center3); + abs3.add_shell(shell3.pure(), shell3.l(), cg3); + ket3.add_center(abs3); + + const auto shell4 = water_aos.shell(4); + const auto center4 = shell4.center().as_point(); + const auto cg4_view = shell4.contracted_gaussian(); + vector_type cs4{cg4_view[0].coefficient(), cg4_view[1].coefficient(), + cg4_view[2].coefficient()}; + vector_type es4{cg4_view[0].exponent(), cg4_view[1].exponent(), + cg4_view[2].exponent()}; + cg_type cg4(cs4.begin(), cs4.end(), es4.begin(), es4.end(), center4); + + abs_type abs4(center4); + abs4.add_shell(shell4.pure(), shell4.l(), cg4); + ket4.add_center(abs4); + return std::make_tuple(bra0, bra0, ket3, ket4); +} + } // namespace integrals::testing \ No newline at end of file diff --git a/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp b/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp index a12747c7..20ebf986 100644 --- a/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp +++ b/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp @@ -17,6 +17,36 @@ #include using namespace integrals::testing; +using pair_type = std::vector; +using pair_vector = std::vector; + +namespace { + +auto corr_bra_pairs() { + pair_vector rv; + rv.push_back({0, 2}); + rv.push_back({1, 1}); + rv.push_back({1, 2}); + rv.push_back({2, 0}); + rv.push_back({2, 1}); + rv.push_back({2, 2}); + return rv; +} + +auto corr_ket_pairs() { + pair_vector rv; + rv.push_back({0, 1}); + rv.push_back({0, 2}); + rv.push_back({1, 0}); + rv.push_back({1, 1}); + rv.push_back({1, 2}); + rv.push_back({2, 0}); + rv.push_back({2, 1}); + rv.push_back({2, 2}); + return rv; +} + +} // namespace using pt = integrals::property_types::PairScreener; @@ -29,9 +59,11 @@ TEST_CASE("Screen Primitive Pairs") { SECTION("(Bra0, Bra1) pairs") { auto pairs = mod.run_as(bra0, bra1, 1E-10); - for(const auto& pair : pairs) { - for(const auto& prim : pair) { std::cout << prim << " "; } - std::cout << std::endl; - } + REQUIRE(pairs == corr_bra_pairs()); + } + + SECTION("(Ket0, Ket1) pairs") { + auto pairs = mod.run_as(ket0, ket1, 1E-10); + REQUIRE(pairs == corr_ket_pairs()); } } \ No newline at end of file From 032709fc5c8310e8502b392b5382db28ffe27b82 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 24 Feb 2026 11:18:53 -0600 Subject: [PATCH 10/14] run precommit/add missing header --- examples/primitive_error_models/main.cpp | 18 +++++++++++++++++- src/integrals/ao_integrals/uq_driver.cpp | 1 + src/integrals/libint/analytic_error.cpp | 18 +++++++++++++++++- src/integrals/libint/primitive_error_model.cpp | 18 +++++++++++++++++- src/integrals/utils/rank2_shell_norm.hpp | 18 +++++++++++++++++- src/integrals/utils/screen_primitive_pairs.cpp | 18 +++++++++++++++++- .../unit/integrals/libint/analytic_error.cpp | 18 +++++++++++++++++- .../integrals/libint/primitive_error_model.cpp | 18 +++++++++++++++++- tests/cxx/unit/integrals/testing/ao_bases.hpp | 18 +++++++++++++++++- tests/cxx/unit/integrals/testing/molecules.hpp | 18 +++++++++++++++++- .../unit/integrals/testing/shell_quartets.hpp | 18 +++++++++++++++++- tests/cxx/unit/integrals/testing/testing.hpp | 18 +++++++++++++++++- 12 files changed, 188 insertions(+), 11 deletions(-) diff --git a/examples/primitive_error_models/main.cpp b/examples/primitive_error_models/main.cpp index bd2aa627..8497d79c 100644 --- a/examples/primitive_error_models/main.cpp +++ b/examples/primitive_error_models/main.cpp @@ -1,3 +1,19 @@ +/* + * 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 #include @@ -71,4 +87,4 @@ int main(int argc, char* argv[]) { std::cout << "Estimated error: " << approx_error << std::endl; return 0; -} \ No newline at end of file +} diff --git a/src/integrals/ao_integrals/uq_driver.cpp b/src/integrals/ao_integrals/uq_driver.cpp index 9d03ba0c..4d464546 100644 --- a/src/integrals/ao_integrals/uq_driver.cpp +++ b/src/integrals/ao_integrals/uq_driver.cpp @@ -16,6 +16,7 @@ #include "ao_integrals.hpp" #include +#include using namespace tensorwrapper; diff --git a/src/integrals/libint/analytic_error.cpp b/src/integrals/libint/analytic_error.cpp index e4b01005..33420b4e 100644 --- a/src/integrals/libint/analytic_error.cpp +++ b/src/integrals/libint/analytic_error.cpp @@ -1,3 +1,19 @@ +/* + * 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 "libint.hpp" #include @@ -37,4 +53,4 @@ MODULE_RUN(AnalyticError) { return pt::wrap_results(rv, error); } -} // namespace integrals::libint \ No newline at end of file +} // namespace integrals::libint diff --git a/src/integrals/libint/primitive_error_model.cpp b/src/integrals/libint/primitive_error_model.cpp index 6cdce089..55749c4f 100644 --- a/src/integrals/libint/primitive_error_model.cpp +++ b/src/integrals/libint/primitive_error_model.cpp @@ -1,3 +1,19 @@ +/* + * 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 "libint.hpp" #include #include @@ -201,4 +217,4 @@ MODULE_RUN(PrimitiveErrorModel) { return pt::wrap_results(result, error); } -} // namespace integrals::libint \ No newline at end of file +} // namespace integrals::libint diff --git a/src/integrals/utils/rank2_shell_norm.hpp b/src/integrals/utils/rank2_shell_norm.hpp index e97fa252..84770dc9 100644 --- a/src/integrals/utils/rank2_shell_norm.hpp +++ b/src/integrals/utils/rank2_shell_norm.hpp @@ -1,3 +1,19 @@ +/* + * 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 @@ -38,4 +54,4 @@ auto rank2_shell_norm(TensorType&& buffer, std::array offsets, return std::sqrt(shell_norm); } -} // namespace integrals::utils \ No newline at end of file +} // namespace integrals::utils diff --git a/src/integrals/utils/screen_primitive_pairs.cpp b/src/integrals/utils/screen_primitive_pairs.cpp index a4cc06c8..51832c9c 100644 --- a/src/integrals/utils/screen_primitive_pairs.cpp +++ b/src/integrals/utils/screen_primitive_pairs.cpp @@ -1,3 +1,19 @@ +/* + * 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 "utils.hpp" #include namespace integrals::utils { @@ -55,4 +71,4 @@ MODULE_RUN(ScreenPrimitivePairs) { return pt::wrap_results(result, rv); } -} // namespace integrals::utils \ No newline at end of file +} // namespace integrals::utils diff --git a/tests/cxx/unit/integrals/libint/analytic_error.cpp b/tests/cxx/unit/integrals/libint/analytic_error.cpp index 071d9889..57c4fc8f 100644 --- a/tests/cxx/unit/integrals/libint/analytic_error.cpp +++ b/tests/cxx/unit/integrals/libint/analytic_error.cpp @@ -1,3 +1,19 @@ +/* + * 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 "../testing/testing.hpp" #include @@ -48,4 +64,4 @@ TEST_CASE("AnalyticError") { auto corr = corr_error(-0.0000000000315610); REQUIRE(approximately_equal(error, corr, 1.0e-12)); } -} \ No newline at end of file +} diff --git a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp index efd91718..cc3efe1f 100644 --- a/tests/cxx/unit/integrals/libint/primitive_error_model.cpp +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -1,3 +1,19 @@ +/* + * 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 "../testing/testing.hpp" #include @@ -96,4 +112,4 @@ TEST_CASE("PrimitiveErrorModel") { auto corr = anal_error.run_as(mnls, tol); REQUIRE(compare_magnitudes(error, corr, tol)); } -} \ No newline at end of file +} diff --git a/tests/cxx/unit/integrals/testing/ao_bases.hpp b/tests/cxx/unit/integrals/testing/ao_bases.hpp index e4a619cd..f423fbb0 100644 --- a/tests/cxx/unit/integrals/testing/ao_bases.hpp +++ b/tests/cxx/unit/integrals/testing/ao_bases.hpp @@ -1,3 +1,19 @@ +/* + * 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 "molecules.hpp" #include @@ -126,4 +142,4 @@ inline simde::type::ao_basis_set water_decontracted_sto3g_basis_set() { return bs; } -} // namespace integrals::testing \ No newline at end of file +} // namespace integrals::testing diff --git a/tests/cxx/unit/integrals/testing/molecules.hpp b/tests/cxx/unit/integrals/testing/molecules.hpp index ac87f77d..5a140db8 100644 --- a/tests/cxx/unit/integrals/testing/molecules.hpp +++ b/tests/cxx/unit/integrals/testing/molecules.hpp @@ -1,3 +1,19 @@ +/* + * 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 @@ -22,4 +38,4 @@ inline simde::type::molecule water_molecule() { return molecule_t{O, H1, H2}; } -} // namespace integrals::testing \ No newline at end of file +} // namespace integrals::testing diff --git a/tests/cxx/unit/integrals/testing/shell_quartets.hpp b/tests/cxx/unit/integrals/testing/shell_quartets.hpp index 48d99085..abf0cbdd 100644 --- a/tests/cxx/unit/integrals/testing/shell_quartets.hpp +++ b/tests/cxx/unit/integrals/testing/shell_quartets.hpp @@ -1,3 +1,19 @@ +/* + * 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 "ao_bases.hpp" #include @@ -86,4 +102,4 @@ inline auto get_h2o_0034_bases() { return std::make_tuple(bra0, bra0, ket3, ket4); } -} // namespace integrals::testing \ No newline at end of file +} // namespace integrals::testing diff --git a/tests/cxx/unit/integrals/testing/testing.hpp b/tests/cxx/unit/integrals/testing/testing.hpp index cb6a7e3b..035b30b1 100644 --- a/tests/cxx/unit/integrals/testing/testing.hpp +++ b/tests/cxx/unit/integrals/testing/testing.hpp @@ -1,3 +1,19 @@ +/* + * 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. + */ + /* Copyright 2024 NWChemEx - Project * Licensed under the Apache License, Version 2.0(the "License"); * you may not use this file except in compliance with the License. @@ -79,4 +95,4 @@ inline auto h2_density() { return density_type(rho, h2_mos()); } -} // namespace integrals::testing \ No newline at end of file +} // namespace integrals::testing From 10bf4a73b6be6d8209b2be5d5eb0da8517547c02 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 24 Feb 2026 11:49:58 -0600 Subject: [PATCH 11/14] ifdef protect sigma code --- src/integrals/ao_integrals/uq_driver.cpp | 6 ++++++ src/integrals/libint/black_box_primitive_estimator.cpp | 8 ++++---- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/integrals/ao_integrals/uq_driver.cpp b/src/integrals/ao_integrals/uq_driver.cpp index 4d464546..43ee1f87 100644 --- a/src/integrals/ao_integrals/uq_driver.cpp +++ b/src/integrals/ao_integrals/uq_driver.cpp @@ -16,7 +16,9 @@ #include "ao_integrals.hpp" #include +#ifdef ENABLE_SIGMA #include +#endif using namespace tensorwrapper; @@ -43,6 +45,7 @@ struct Kernel { if constexpr(types::is_uncertain_v) { throw std::runtime_error("Did not expect an uncertain type"); } else { +#ifdef ENABLE_SIGMA using uq_type = sigma::Uncertain; auto rv_buffer = buffer::make_contiguous(m_shape); auto rv_data = buffer::get_raw_data(rv_buffer); @@ -52,6 +55,9 @@ struct Kernel { rv_data[i] = uq_type(elem, elem_error); } rv = tensorwrapper::Tensor(m_shape, std::move(rv_buffer)); +#else + throw std::runtime_error("Sigma support not enabled!"); +#endif } return rv; diff --git a/src/integrals/libint/black_box_primitive_estimator.cpp b/src/integrals/libint/black_box_primitive_estimator.cpp index 2c76c25d..6720c8a0 100644 --- a/src/integrals/libint/black_box_primitive_estimator.cpp +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -23,14 +23,14 @@ namespace integrals::libint { namespace { const auto desc = R"( -Libint Black Box Primitive Pair Estimator +Libint Black Box Primitive Pair Estimator ========================================= This module computes the matrix : math :`K_{ij}` where: -.. math:: +.. math:: - K_{ij} = c_i c_j \exp\left(-\frac{\zeta_i \zeta_j}{\zeta_i + \zeta_j} + K_{ij} = c_i c_j \exp\left(-\frac{\zeta_i \zeta_j}{\zeta_i + \zeta_j} |\mathbf{R}_i - \mathbf{R}_j|^2\right) This is how Libint2 estimates the contribution of a pair of primitives to @@ -154,4 +154,4 @@ MODULE_RUN(BlackBoxPrimitiveEstimator) { return pt::wrap_results(result, rv); } -} // namespace integrals::libint \ No newline at end of file +} // namespace integrals::libint From 4fe81d9b6d3e863a0b204f452eaa06185a258996 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 24 Feb 2026 13:11:10 -0600 Subject: [PATCH 12/14] blank lines for precommit --- src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp | 2 +- .../cxx/unit/integrals/libint/black_box_primitive_estimator.cpp | 2 +- .../integrals/libint/cauchy_schwarz_primitive_estimator.cpp | 2 +- tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp | 2 +- tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp index c6665ca0..2a6718ea 100644 --- a/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp +++ b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -146,4 +146,4 @@ MODULE_RUN(CauchySchwarzPrimitiveEstimator) { return pt::wrap_results(result, rv); } -} // namespace integrals::libint \ No newline at end of file +} // namespace integrals::libint diff --git a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp index 3eb07127..43ed6671 100644 --- a/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp +++ b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp @@ -116,4 +116,4 @@ TEST_CASE("BlackBoxPrimitiveEstimator") { REQUIRE(approximately_equal(K23, corr, 1E-8)); } } -} \ No newline at end of file +} diff --git a/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp index 9f297a5b..072e53e8 100644 --- a/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp +++ b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -105,4 +105,4 @@ TEST_CASE("CauchySchwarzPrimitiveEstimator") { auto corr = Q_ab_pppp_corr(); REQUIRE(approximately_equal(Q_ab, corr, 1E-8)); } -} \ No newline at end of file +} diff --git a/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp b/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp index 07657803..445eff80 100644 --- a/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp +++ b/tests/cxx/unit/integrals/utils/rank2_shell_norm.cpp @@ -96,4 +96,4 @@ TEST_CASE("Rank2 Shell Norm") { auto norm = rank2_shell_norm(buffer, offsets, naos); REQUIRE(norm == Approx(0.0035472410738408154).epsilon(1.0e-8)); } -} \ No newline at end of file +} diff --git a/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp b/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp index 20ebf986..422f8aa0 100644 --- a/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp +++ b/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp @@ -66,4 +66,4 @@ TEST_CASE("Screen Primitive Pairs") { auto pairs = mod.run_as(ket0, ket1, 1E-10); REQUIRE(pairs == corr_ket_pairs()); } -} \ No newline at end of file +} From a5609d21eaa753c167c5da16ab75bc92b74aac8a Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 24 Feb 2026 14:35:11 -0600 Subject: [PATCH 13/14] address comments --- src/integrals/libint/black_box_primitive_estimator.cpp | 8 -------- src/integrals/utils/utils.hpp | 1 - 2 files changed, 9 deletions(-) diff --git a/src/integrals/libint/black_box_primitive_estimator.cpp b/src/integrals/libint/black_box_primitive_estimator.cpp index 6720c8a0..53717901 100644 --- a/src/integrals/libint/black_box_primitive_estimator.cpp +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -59,14 +59,6 @@ auto compute_k(T zeta_i, T zeta_j, T coeff_i, T coeff_j, T dr2) { return coeff_i * coeff_j * std::exp(ratio * dr2); } -template -auto compute_ln_k(T zeta_i, T zeta_j, T coeff_i, T coeff_j, T dr2) { - const auto num = -zeta_i * zeta_j; - const auto denom = zeta_i + zeta_j; - const auto ratio = num / denom; - return std::log(coeff_i) + std::log(coeff_j) + ratio * dr2; -} - } // namespace using pt = integrals::property_types::PrimitivePairEstimator; diff --git a/src/integrals/utils/utils.hpp b/src/integrals/utils/utils.hpp index 72e3cc36..f2af1a56 100644 --- a/src/integrals/utils/utils.hpp +++ b/src/integrals/utils/utils.hpp @@ -40,7 +40,6 @@ inline void set_defaults(pluginplay::ModuleManager& mm) { inline void load_modules(pluginplay::ModuleManager& mm) { mm.add_module("Decontract Basis Set"); mm.add_module("Screen Primitive Pairs"); - set_defaults(mm); } } // end namespace integrals::utils From 9f1f55d05bf0b44825c79542a37ff0afe95ad12e Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Tue, 24 Feb 2026 15:03:36 -0600 Subject: [PATCH 14/14] call utils::set_defaults --- src/integrals/integrals_mm.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/integrals/integrals_mm.cpp b/src/integrals/integrals_mm.cpp index 7021528c..9ec4deb1 100644 --- a/src/integrals/integrals_mm.cpp +++ b/src/integrals/integrals_mm.cpp @@ -30,6 +30,7 @@ namespace integrals { void set_defaults(pluginplay::ModuleManager& mm) { libint::set_defaults(mm); ao_integrals::set_defaults(mm); + utils::set_defaults(mm); mm.change_submod("AO integral driver", "Kinetic", "Kinetic"); mm.change_submod("AO integral driver", "Electron-Nuclear attraction", "Nuclear");