diff --git a/.github/enable_sigma.cmake b/.github/enable_sigma.cmake new file mode 100644 index 00000000..2fbe8f6a --- /dev/null +++ b/.github/enable_sigma.cmake @@ -0,0 +1 @@ +set(ENABLE_SIGMA ON) \ No newline at end of file diff --git a/.github/workflows/pull_request.yaml b/.github/workflows/pull_request.yaml index f08bc7f1..f5da17e7 100644 --- a/.github/workflows/pull_request.yaml +++ b/.github/workflows/pull_request.yaml @@ -29,3 +29,13 @@ jobs: compilers: '["gcc-11", "clang-14"]' doc_target: 'integrals_cxx_api' secrets: inherit + Test-Enable-Sigma: + uses: NWChemEx/.github/.github/workflows/common_pull_request.yaml@master + with: + config_file: '' + source_dir: '' + format_python: false + compilers: '["gcc-11", "clang-14"]' + repo_toolchain: ".github/enable_sigma.cmake" + doc_target: '' + secrets: inherit diff --git a/CMakeLists.txt b/CMakeLists.txt index c55816cf..805b326f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ nwx_cxx_api_docs("${INTEGRALS_SOURCE_DIR}" "${INTEGRALS_INCLUDE_DIR}") cmaize_option_list( BUILD_TESTING OFF "Should we build the tests?" BUILD_PYBIND11_PYBINDINGS ON "Build pybind11 python3 bindings?" + ENABLE_SIGMA OFF "Should we enable Sigma for uncertainty tracking?" ) # Can't build from github due to Libint setup. @@ -56,6 +57,7 @@ cmaize_find_or_build_dependency( FIND_TARGET nwx::simde CMAKE_ARGS BUILD_TESTING=OFF BUILD_PYBIND11_PYBINDINGS=${BUILD_PYBIND11_PYBINDINGS} + BUILD_SIGMA=${BUILD_SIGMA} ) cmaize_add_library( diff --git a/src/integrals/ao_integrals/ao_integrals.cpp b/src/integrals/ao_integrals/ao_integrals.cpp index b477759c..8d66aaeb 100644 --- a/src/integrals/ao_integrals/ao_integrals.cpp +++ b/src/integrals/ao_integrals/ao_integrals.cpp @@ -13,99 +13,52 @@ * See the License for the specific language governing permissions and * limitations under the License. */ - +#include "../uncertain_types.hpp" #include "ao_integrals.hpp" #include "detail_/get_basis_sets.hpp" #include "detail_/libint_op.hpp" #include "detail_/make_engine.hpp" #include "detail_/make_libint_basis_set.hpp" #include "detail_/shells2ord.hpp" +#include "lib_int_visitor.hpp" #include namespace integrals::ao_integrals { - -class LibIntVisitor : public chemist::qm_operator::OperatorVisitor { -public: - using s_e_type = simde::type::s_e_type; - using t_e_type = simde::type::t_e_type; - using v_ee_type = simde::type::v_ee_type; - using v_en_type = simde::type::v_en_type; - - LibIntVisitor(const std::vector& bases, double thresh, - std::size_t deriv = 0) : - m_bases(bases), m_thresh(thresh), m_deriv(deriv){}; - - void run(const s_e_type& S_e) { - m_engine = detail_::make_engine(m_bases, S_e, m_thresh, m_deriv); - } - - void run(const t_e_type& T_e) { - m_engine = detail_::make_engine(m_bases, T_e, m_thresh, m_deriv); - } - - void run(const v_en_type& V_en) { - m_engine = detail_::make_engine(m_bases, V_en, m_thresh, m_deriv); - } - - void run(const v_ee_type& V_ee) { - m_engine = detail_::make_engine(m_bases, V_ee, m_thresh, m_deriv); +namespace { + +template +auto build_eigen_buffer(const std::vector& basis_sets, + double thresh) { + FloatType initial_value; + if constexpr(std::is_same_v) { + initial_value = 0.0; + } else { // Presumably sigma::UDouble + initial_value = FloatType(0.0, thresh); } - - libint2::Engine& engine() { return m_engine; } - -private: - const std::vector& m_bases; - double m_thresh; - std::size_t m_deriv; - libint2::Engine m_engine; -}; - -template -TEMPLATED_MODULE_CTOR(AOIntegral, BraKetType) { - using my_pt = simde::EvaluateBraKet; - satisfies_property_type(); - description("Computes integrals with Libint"); - - add_input("Threshold") - .set_default(1.0E-16) - .set_description( - "The target precision with which the integrals will be computed"); -} - -template -TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) { - using my_pt = simde::EvaluateBraKet; - - const auto& [braket] = my_pt::unwrap_inputs(inputs); - auto thresh = inputs.at("Threshold").value(); - auto bra = braket.bra(); - auto ket = braket.ket(); - auto& op = braket.op(); - - // Gather information from Bra, Ket, and Op - auto basis_sets = detail_::get_basis_sets(bra, ket); - constexpr int N = detail_::get_n(bra, ket); - - // Dimensional information - std::vector dims_shells(N); Eigen::array dims_bfs; - for(auto i = 0; i < N; ++i) { - dims_shells[i] = basis_sets[i].size(); - dims_bfs[i] = basis_sets[i].nbf(); - } + for(decltype(N) i = 0; i < N; ++i) dims_bfs[i] = basis_sets[i].nbf(); - // Build tensor inputs - using tensor_t = simde::type::tensor; using shape_t = tensorwrapper::shape::Smooth; using layout_t = tensorwrapper::layout::Physical; - using buffer_t = tensorwrapper::buffer::Eigen; + using buffer_t = tensorwrapper::buffer::Eigen; using data_t = typename buffer_t::data_type; shape_t s{dims_bfs.begin(), dims_bfs.end()}; layout_t l(s); data_t d(dims_bfs); buffer_t b{d, l}; - b.value().setZero(); + b.value().setConstant(initial_value); + return b; +} + +template +auto fill_tensor(const std::vector& basis_sets, + const chemist::qm_operator::OperatorBase& op, double thresh) { + // Dimensional information + std::vector dims_shells(N); + for(decltype(N) i = 0; i < N; ++i) dims_shells[i] = basis_sets[i].size(); + + auto b = build_eigen_buffer(basis_sets, thresh); // Make libint engine LibIntVisitor visitor(basis_sets, thresh); @@ -124,13 +77,13 @@ TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) { auto ord = detail_::shells2ord(basis_sets, shells); auto n_ord = ord.size(); for(decltype(n_ord) i_ord = 0; i_ord < n_ord; ++i_ord) { - b.value().data()[ord[i_ord]] = vals[i_ord]; + b.value().data()[ord[i_ord]] += vals[i_ord]; } } // Increment index shells[N - 1] += 1; - for(auto i = 1; i < N; ++i) { + for(decltype(N) i = 1; i < N; ++i) { if(shells[N - i] >= dims_shells[N - i]) { // Reset this dimension and increment the next one // shells[0] accumulates until we reach the end @@ -140,7 +93,51 @@ TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) { } } - tensor_t t(s, b); + return simde::type::tensor(b.layout().shape().clone(), b); +} + +} // namespace + +template +TEMPLATED_MODULE_CTOR(AOIntegral, BraKetType) { + using my_pt = simde::EvaluateBraKet; + satisfies_property_type(); + description("Driver for computing integrals with Libint"); + + add_input("Threshold") + .set_default(1.0E-16) + .set_description( + "The target precision with which the integrals will be computed"); + + add_input("With UQ?").set_default(false); +} + +template +TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) { + using my_pt = simde::EvaluateBraKet; + + const auto& [braket] = my_pt::unwrap_inputs(inputs); + auto thresh = inputs.at("Threshold").value(); + auto with_uq = inputs.at("With UQ?").value(); + auto bra = braket.bra(); + auto ket = braket.ket(); + auto& op = braket.op(); + + // Gather information from Bra, Ket, and Op + auto basis_sets = detail_::get_basis_sets(bra, ket); + constexpr int N = detail_::get_n(bra, ket); + + simde::type::tensor t; + if(with_uq) { + if constexpr(integrals::type::has_sigma()) { + t = fill_tensor(basis_sets, op, thresh); + } else { + throw std::runtime_error("Sigma support not enabled!"); + } + } else { + t = fill_tensor(basis_sets, op, thresh); + } + auto rv = results(); return my_pt::wrap_results(rv, t); } diff --git a/src/integrals/ao_integrals/lib_int_visitor.hpp b/src/integrals/ao_integrals/lib_int_visitor.hpp new file mode 100644 index 00000000..1ddb2564 --- /dev/null +++ b/src/integrals/ao_integrals/lib_int_visitor.hpp @@ -0,0 +1,57 @@ +/* + * 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 + +namespace integrals::ao_integrals { + +class LibIntVisitor : public chemist::qm_operator::OperatorVisitor { +public: + using s_e_type = simde::type::s_e_type; + using t_e_type = simde::type::t_e_type; + using v_ee_type = simde::type::v_ee_type; + using v_en_type = simde::type::v_en_type; + + LibIntVisitor(const std::vector& bases, double thresh, + std::size_t deriv = 0) : + m_bases(bases), m_thresh(thresh), m_deriv(deriv){}; + + void run(const s_e_type& S_e) { + m_engine = detail_::make_engine(m_bases, S_e, m_thresh, m_deriv); + } + + void run(const t_e_type& T_e) { + m_engine = detail_::make_engine(m_bases, T_e, m_thresh, m_deriv); + } + + void run(const v_en_type& V_en) { + m_engine = detail_::make_engine(m_bases, V_en, m_thresh, m_deriv); + } + + void run(const v_ee_type& V_ee) { + m_engine = detail_::make_engine(m_bases, V_ee, m_thresh, m_deriv); + } + + libint2::Engine& engine() { return m_engine; } + +private: + const std::vector& m_bases; + double m_thresh; + std::size_t m_deriv; + libint2::Engine m_engine; +}; + +} // namespace integrals::ao_integrals \ No newline at end of file diff --git a/src/integrals/uncertain_types.hpp b/src/integrals/uncertain_types.hpp new file mode 100644 index 00000000..f7160372 --- /dev/null +++ b/src/integrals/uncertain_types.hpp @@ -0,0 +1,38 @@ +/* + * Copyright 2025 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 +#ifdef ENABLE_SIGMA +#include +#endif + +namespace integrals::type { +#ifdef ENABLE_SIGMA + +static constexpr bool has_sigma() { return true; } + +using uncertain_float = sigma::UFloat; +using uncertain_double = sigma::UDouble; +#else + +static constexpr bool has_sigma() { return false; } + +using uncertain_float = float; +using uncertain_double = double; +#endif + +} // namespace integrals::type \ No newline at end of file diff --git a/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp b/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp index d11b1a15..7a51fccc 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp @@ -21,22 +21,22 @@ namespace test { -template +template auto eigen_buffer(const tensorwrapper::buffer::BufferBase& buffer) { - return static_cast&>(buffer); + return static_cast&>( + buffer); } -template -auto trace(const TensorType& t) { - Eigen::Tensor trace = t.value().trace(); +template +auto trace(const tensorwrapper::buffer::Eigen& t) { + Eigen::Tensor trace = t.value().trace(); return trace.coeff(); } -template -auto norm(const TensorType& t) { - Eigen::Tensor norm = +template +auto norm(const tensorwrapper::buffer::Eigen& t) { + Eigen::Tensor norm = t.value().square().sum().sqrt(); return norm.coeff(); } - } // namespace test \ No newline at end of file diff --git a/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp b/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp index 1b63467f..f20d312d 100644 --- a/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp +++ b/tests/cxx/unit/integrals/ao_integrals/test_arbitrary_operator.cpp @@ -14,8 +14,30 @@ * limitations under the License. */ +#include "integrals/uncertain_types.hpp" #include "test_ao_integrals.hpp" +using udouble = integrals::type::uncertain_double; +constexpr bool has_sigma = integrals::type::has_sigma(); + +template +auto unwrap_mean(InputType&& uq) { + if constexpr(has_sigma) { + return uq.mean(); + } else { + return uq; + } +} + +template +auto unwrap_sd(InputType&& uq) { + if constexpr(has_sigma) { + return uq.sd(); + } else { + return 0.0; + } +} + TEST_CASE("OperatorBase") { using aos_t = simde::type::aos; using aos_squared_t = simde::type::aos_squared; @@ -48,16 +70,37 @@ TEST_CASE("OperatorBase") { braket_t braket(aos, op_base, aos); // Call module - auto T = mm.at("Evaluate 2-Index BraKet").run_as(braket); - - // Check output - auto t = test::eigen_buffer<2>(T.buffer()); - REQUIRE(test::trace(t) == - Catch::Approx(124.7011973877891364).margin(1.0e-16)); - REQUIRE(test::norm(t) == - Catch::Approx(90.2562579028763707).margin(1.0e-16)); + auto& mod = mm.at("Evaluate 2-Index BraKet"); + + SECTION("No UQ") { + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<2>(T.buffer()); + REQUIRE(test::trace(t) == + Catch::Approx(124.7011973877891364).margin(1.0e-16)); + REQUIRE(test::norm(t) == + Catch::Approx(90.2562579028763707).margin(1.0e-16)); + } + + SECTION("With UQ") { + if constexpr(has_sigma) { + mod.change_input("With UQ?", true); + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<2, udouble>(T.buffer()); + REQUIRE(unwrap_mean(test::trace(t)) == + Catch::Approx(124.7011973877891364).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::trace(t)) == + Catch::Approx(7e-16).margin(1.0e-16)); + REQUIRE(unwrap_mean(test::norm(t)) == + Catch::Approx(90.2562579028763707).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::norm(t)) == + Catch::Approx(3e-16).margin(1.0e-16)); + } + } } - SECTION("3-Index") { using braket_t = simde::type::braket; using test_pt = simde::EvaluateBraKet; @@ -65,15 +108,38 @@ TEST_CASE("OperatorBase") { // Make BraKet Input braket_t braket(aos, op_base, aos_squared); - // Call module - auto T = mm.at("Evaluate 3-Index BraKet").run_as(braket); - - // Check output - auto t = test::eigen_buffer<3>(T.buffer()); - REQUIRE(test::trace(t) == - Catch::Approx(16.8245948391706577).margin(1.0e-16)); - REQUIRE(test::norm(t) == - Catch::Approx(20.6560572032543597).margin(1.0e-16)); + auto& mod = mm.at("Evaluate 3-Index BraKet"); + + SECTION("No UQ") { + // Call module + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<3>(T.buffer()); + REQUIRE(test::trace(t) == + Catch::Approx(16.8245948391706577).margin(1.0e-16)); + REQUIRE(test::norm(t) == + Catch::Approx(20.6560572032543597).margin(1.0e-16)); + } + + SECTION("With UQ") { + if constexpr(has_sigma) { + mod.change_input("With UQ?", true); + // Call module + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<3, udouble>(T.buffer()); + REQUIRE(unwrap_mean(test::trace(t)) == + Catch::Approx(16.8245948391706577).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::trace(t)) == + Catch::Approx(7e-16).margin(1.0e-16)); + REQUIRE(unwrap_mean(test::norm(t)) == + Catch::Approx(20.6560572032543597).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::norm(t)) == + Catch::Approx(7e-16).margin(1.0e-16)); + } + } } SECTION("4-Index") { @@ -84,14 +150,37 @@ TEST_CASE("OperatorBase") { // Make BraKet Input braket_t braket(aos_squared, op_base, aos_squared); - // Call module - auto T = mm.at("Evaluate 4-Index BraKet").run_as(braket); - - // Check output - auto t = test::eigen_buffer<4>(T.buffer()); - REQUIRE(test::trace(t) == - Catch::Approx(9.7919608941952063).margin(1.0e-16)); - REQUIRE(test::norm(t) == - Catch::Approx(7.7796143419802553).margin(1.0e-16)); + auto& mod = mm.at("Evaluate 4-Index BraKet"); + + SECTION("No UQ") { + // Call module + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<4>(T.buffer()); + REQUIRE(test::trace(t) == + Catch::Approx(9.7919608941952063).margin(1.0e-16)); + REQUIRE(test::norm(t) == + Catch::Approx(7.7796143419802553).margin(1.0e-16)); + } + + SECTION("With UQ") { + if constexpr(has_sigma) { + // Call module + mod.change_input("With UQ?", true); + auto T = mod.run_as(braket); + + // Check output + auto t = test::eigen_buffer<4, udouble>(T.buffer()); + REQUIRE(unwrap_mean(test::trace(t)) == + Catch::Approx(9.7919608941952063).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::trace(t)) == + Catch::Approx(7e-16).margin(1.0e-16)); + REQUIRE(unwrap_mean(test::norm(t)) == + Catch::Approx(7.7796143419802553).margin(1.0e-16)); + REQUIRE(unwrap_sd(test::norm(t)) == + Catch::Approx(11e-16).margin(1.0e-16)); + } + } } }