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..8497d79c --- /dev/null +++ b/examples/primitive_error_models/main.cpp @@ -0,0 +1,90 @@ +/* + * 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 + +/* 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; +} 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..3998de81 100644 --- a/include/integrals/property_types.hpp +++ b/include/integrals/property_types.hpp @@ -30,6 +30,64 @@ */ 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; +} + +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); + +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..6fe9fb10 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", "Primitive Error Model"); } 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..43ee1f87 100644 --- a/src/integrals/ao_integrals/uq_driver.cpp +++ b/src/integrals/ao_integrals/uq_driver.cpp @@ -15,6 +15,10 @@ */ #include "ao_integrals.hpp" +#include +#ifdef ENABLE_SIGMA +#include +#endif using namespace tensorwrapper; @@ -37,19 +41,25 @@ 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 { +#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); 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"); +#else + throw std::runtime_error("Sigma support not enabled!"); +#endif } + return rv; } shape_type m_shape; @@ -63,35 +73,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..9ec4deb1 100644 --- a/src/integrals/integrals_mm.cpp +++ b/src/integrals/integrals_mm.cpp @@ -28,6 +28,9 @@ namespace integrals { * @throw none No throw guarantee */ 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"); @@ -41,8 +44,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..33420b4e --- /dev/null +++ b/src/integrals/libint/analytic_error.cpp @@ -0,0 +1,56 @@ +/* + * 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 + +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 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..53717901 --- /dev/null +++ b/src/integrals/libint/black_box_primitive_estimator.cpp @@ -0,0 +1,149 @@ +/* + * 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 "detail_/make_libint_basis_set.hpp" +#include "libint.hpp" +#include +#include + +namespace integrals::libint { +namespace { + +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. + +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[0] - b[0]; + auto dy = a[1] - b[1]; + auto dz = a[2] - b[2]; + 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) { + satisfies_property_type(); + description(desc); + // TODO: Add citation for Chemist paper +} + +MODULE_RUN(BlackBoxPrimitiveEstimator) { + const auto&& [bra, ket] = pt::unwrap_inputs(inputs); + + 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 + + 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 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 + + // 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 now use the libint basis sets because they're properly normalized + // TODO: Our basis really needs to handle normalization better... + 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); // No general contraction support + 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(); + + 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; + abs_prim[0] = offsets[0] + prims[0]; + + offsets[1] = 0; + for(shells[1] = 0; shells[1] < n_shells[1]; ++shells[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.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] = 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_elem(abs_prim, k01); + } // loop over ket primitives + + offsets[1] += n_prims_ket_shell; + } // loop over ket shells + + // We ran over all ket primitives, so counter should be done too + assert(offsets[1] == n_prims[1]); + } // loop over bra primitives + + offsets[0] += n_prims_bra_shell; + } // loop over bra shells + + // We ran over all bra primitives, so counter should be done too + assert(offsets[0] == n_prims[0]); + + simde::type::tensor rv(shape, std::move(buffer)); + + auto result = results(); + return pt::wrap_results(result, rv); +} + +} // namespace integrals::libint 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..2a6718ea --- /dev/null +++ b/src/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -0,0 +1,149 @@ +/* + * 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 "detail_/make_libint_basis_set.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_prims = bra_prims.n_primitives(); + const auto n_ket_prims = ket_prims.n_primitives(); + + // Should always be true, but we check for sanity + assert(n_bra_prims == bra_basis.n_primitives()); + assert(n_ket_prims == ket_basis.n_primitives()); + + // 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 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); +} + +} // namespace integrals::libint diff --git a/src/integrals/libint/libint.cpp b/src/integrals/libint/libint.cpp index d391c2a0..57fadc32 100644 --- a/src/integrals/libint/libint.cpp +++ b/src/integrals/libint/libint.cpp @@ -192,7 +192,15 @@ 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", + "Black Box Primitive Pair Estimator", + "Black Box Primitive Pair Estimator"); + mm.change_submod("CauchySchwarz Estimator", "Decontract Basis Set", + "Decontract Basis Set"); + mm.copy_module("ERI4", "Benchmark ERI4"); + mm.change_input("Benchmark ERI4", "Threshold", 1.0E-16); + mm.change_submod("CauchySchwarz Estimator", "ERI4", "Benchmark ERI4"); + mm.change_submod("Analytic Error", "ERI4s", "Benchmark ERI4"); } #define LOAD_LIBINT(bra, op, ket, key) mm.add_module(key) @@ -208,7 +216,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..55749c4f --- /dev/null +++ b/src/integrals/libint/primitive_error_model.cpp @@ -0,0 +1,220 @@ +/* + * 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 = "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&& 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(); + 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 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; + + 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 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 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 += tol; + } + } // 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("Black Box 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 K_ab and K_cd tensors (used to determine which primitives are + // neglected) + 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& 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(); + 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, 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] + + 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); +} + +} // namespace integrals::libint diff --git a/src/integrals/utils/rank2_shell_norm.hpp b/src/integrals/utils/rank2_shell_norm.hpp new file mode 100644 index 00000000..84770dc9 --- /dev/null +++ b/src/integrals/utils/rank2_shell_norm.hpp @@ -0,0 +1,57 @@ +/* + * 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 +#include +#include + +namespace integrals::utils { + +template +auto rank2_shell_norm(TensorType&& buffer, std::array offsets, + std::array naos) { + using float_type = double; + + 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 diff --git a/src/integrals/utils/screen_primitive_pairs.cpp b/src/integrals/utils/screen_primitive_pairs.cpp new file mode 100644 index 00000000..51832c9c --- /dev/null +++ b/src/integrals/utils/screen_primitive_pairs.cpp @@ -0,0 +1,74 @@ +/* + * 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 { + +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 diff --git a/src/integrals/utils/utils.hpp b/src/integrals/utils/utils.hpp index 5109d29e..f2af1a56 100644 --- a/src/integrals/utils/utils.hpp +++ b/src/integrals/utils/utils.hpp @@ -30,12 +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"); - set_defaults(mm); + mm.add_module("Screen Primitive Pairs"); } } // end namespace integrals::utils 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..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,21 +52,20 @@ TEST_CASE("AOIntegralsDriver") { chemist::braket::BraKet; - pluginplay::ModuleManager mm; - integrals::load_modules(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 6ceea87d..031c59ff 100644 --- a/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/coulomb_metric.cpp @@ -14,18 +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); + 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); @@ -39,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 e3045bf4..7a36542e 100644 --- a/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/df_integral.cpp @@ -14,18 +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); + 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); @@ -40,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 95777064..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,24 +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); + 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); @@ -40,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 c455391b..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,24 +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); + 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); @@ -39,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 f2341506..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,24 +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); + 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); @@ -40,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 97da531e..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,24 +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); + 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); @@ -39,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 055f1280..e85ab324 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; @@ -53,13 +55,15 @@ 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.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 - 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); @@ -71,10 +75,15 @@ TEST_CASE("UQ Driver") { // Make BraKet Input chemist::braket::BraKet braket(aos_squared, op, aos_squared); - // Call module + // Call modules 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_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/analytic_error.cpp b/tests/cxx/unit/integrals/libint/analytic_error.cpp new file mode 100644 index 00000000..57c4fc8f --- /dev/null +++ b/tests/cxx/unit/integrals/libint/analytic_error.cpp @@ -0,0 +1,67 @@ +/* + * 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 eri4_pt = simde::ERI4; +using pt = integrals::property_types::Uncertainty; + +using namespace integrals::testing; +using tensorwrapper::operations::approximately_equal; +namespace { + +template +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); + 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{}; + 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); + 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)); + } +} 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..43ed6671 --- /dev/null +++ b/tests/cxx/unit/integrals/libint/black_box_primitive_estimator.cpp @@ -0,0 +1,119 @@ +/* + * 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 tensorwrapper::operations::approximately_equal; + +using pt = integrals::property_types::PrimitivePairEstimator; + +namespace { + +/* 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. + */ + +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); +} + +/* 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"); + using float_type = double; + + 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("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)); + } + } +} 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..072e53e8 --- /dev/null +++ b/tests/cxx/unit/integrals/libint/cauchy_schwarz_primitive_estimator.cpp @@ -0,0 +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; +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)); + } +} 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 new file mode 100644 index 00000000..cc3efe1f --- /dev/null +++ b/tests/cxx/unit/integrals/libint/primitive_error_model.cpp @@ -0,0 +1,115 @@ +/* + * 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 eri4_pt = simde::ERI4; +using pt = integrals::property_types::Uncertainty; + +using namespace integrals::testing; + +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 + +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{}; + + 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)); + } +} 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_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 63% rename from tests/cxx/unit/integrals/testing.hpp rename to tests/cxx/unit/integrals/testing/ao_bases.hpp index 812b05c3..f423fbb0 100644 --- a/tests/cxx/unit/integrals/testing.hpp +++ b/tests/cxx/unit/integrals/testing/ao_bases.hpp @@ -1,5 +1,5 @@ /* - * Copyright 2024 NWChemEx-Project + * 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. @@ -15,60 +15,35 @@ */ #pragma once -#include - -// For common inputs +#include "molecules.hpp" #include -// Common Catch2 includes -#include -#include - -#include - -namespace test { - -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() { @@ -167,53 +142,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; - - 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 diff --git a/tests/cxx/unit/integrals/testing/molecules.hpp b/tests/cxx/unit/integrals/testing/molecules.hpp new file mode 100644 index 00000000..5a140db8 --- /dev/null +++ b/tests/cxx/unit/integrals/testing/molecules.hpp @@ -0,0 +1,41 @@ +/* + * 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 + +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 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..abf0cbdd --- /dev/null +++ b/tests/cxx/unit/integrals/testing/shell_quartets.hpp @@ -0,0 +1,105 @@ +/* + * 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 + +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); +} + +/* 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 diff --git a/tests/cxx/unit/integrals/testing/testing.hpp b/tests/cxx/unit/integrals/testing/testing.hpp new file mode 100644 index 00000000..035b30b1 --- /dev/null +++ b/tests/cxx/unit/integrals/testing/testing.hpp @@ -0,0 +1,98 @@ +/* + * 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. + * 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 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..445eff80 --- /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 "../testing/testing.hpp" +#include + +using namespace integrals::testing; +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)); + } +} 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..422f8aa0 --- /dev/null +++ b/tests/cxx/unit/integrals/utils/screen_primitive_pairs.cpp @@ -0,0 +1,69 @@ +/* + * 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 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; + +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); + REQUIRE(pairs == corr_bra_pairs()); + } + + SECTION("(Ket0, Ket1) pairs") { + auto pairs = mod.run_as(ket0, ket1, 1E-10); + REQUIRE(pairs == corr_ket_pairs()); + } +} 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); }