Adds Sigma

.github/enable_sigma.cmake

@@ -0,0 +1 @@
+set(ENABLE_SIGMA ON)

.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

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(

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 <type_traits>
 
 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<libint2::BasisSet>& 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<typename FloatType, unsigned int N>
+auto build_eigen_buffer(const std::vector<libint2::BasisSet>& basis_sets,
+                        double thresh) {
+    FloatType initial_value;
+    if constexpr(std::is_same_v<FloatType, double>) {
+        initial_value = 0.0;
+    } else { // Presumably sigma::UDouble
+        initial_value = FloatType(0.0, thresh);
     }
-
-    libint2::Engine& engine() { return m_engine; }
-
-private:
-    const std::vector<libint2::BasisSet>& m_bases;
-    double m_thresh;
-    std::size_t m_deriv;
-    libint2::Engine m_engine;
-};
-
-template<typename BraKetType>
-TEMPLATED_MODULE_CTOR(AOIntegral, BraKetType) {
-    using my_pt = simde::EvaluateBraKet<BraKetType>;
-    satisfies_property_type<my_pt>();
-    description("Computes integrals with Libint");
-
-    add_input<double>("Threshold")
-      .set_default(1.0E-16)
-      .set_description(
-        "The target precision with which the integrals will be computed");
-}
-
-template<typename BraKetType>
-TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) {
-    using my_pt = simde::EvaluateBraKet<BraKetType>;
-
-    const auto& [braket] = my_pt::unwrap_inputs(inputs);
-    auto thresh          = inputs.at("Threshold").value<double>();
-    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<std::size_t> dims_shells(N);
     Eigen::array<Eigen::Index, N> 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<double, N>;
+    using buffer_t = tensorwrapper::buffer::Eigen<FloatType, N>;
     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<typename FloatType, unsigned int N>
+auto fill_tensor(const std::vector<libint2::BasisSet>& basis_sets,
+                 const chemist::qm_operator::OperatorBase& op, double thresh) {
+    // Dimensional information
+    std::vector<std::size_t> dims_shells(N);
+    for(decltype(N) i = 0; i < N; ++i) dims_shells[i] = basis_sets[i].size();
+
+    auto b = build_eigen_buffer<FloatType, N>(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<typename BraKetType>
+TEMPLATED_MODULE_CTOR(AOIntegral, BraKetType) {
+    using my_pt = simde::EvaluateBraKet<BraKetType>;
+    satisfies_property_type<my_pt>();
+    description("Driver for computing integrals with Libint");
+
+    add_input<double>("Threshold")
+      .set_default(1.0E-16)
+      .set_description(
+        "The target precision with which the integrals will be computed");
+
+    add_input<bool>("With UQ?").set_default(false);
+}
+
+template<typename BraKetType>
+TEMPLATED_MODULE_RUN(AOIntegral, BraKetType) {
+    using my_pt = simde::EvaluateBraKet<BraKetType>;
+
+    const auto& [braket] = my_pt::unwrap_inputs(inputs);
+    auto thresh          = inputs.at("Threshold").value<double>();
+    auto with_uq         = inputs.at("With UQ?").value<bool>();
+    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<type::uncertain_double, N>(basis_sets, op, thresh);
+        } else {
+            throw std::runtime_error("Sigma support not enabled!");
+        }
+    } else {
+        t = fill_tensor<double, N>(basis_sets, op, thresh);
+    }
+
     auto rv = results();
     return my_pt::wrap_results(rv, t);
 }

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 <simde/simde.hpp>
+
+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<libint2::BasisSet>& 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<libint2::BasisSet>& m_bases;
+    double m_thresh;
+    std::size_t m_deriv;
+    libint2::Engine m_engine;
+};
+
+} // namespace integrals::ao_integrals

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 <stdexcept>
+#ifdef ENABLE_SIGMA
+#include <sigma/sigma.hpp>
+#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

tests/cxx/unit/integrals/ao_integrals/test_ao_integrals.hpp

@@ -21,22 +21,22 @@
 
 namespace test {
 
-template<std::size_t N>
+template<std::size_t N, typename FloatType = double>
 auto eigen_buffer(const tensorwrapper::buffer::BufferBase& buffer) {
-    return static_cast<const tensorwrapper::buffer::Eigen<double, N>&>(buffer);
+    return static_cast<const tensorwrapper::buffer::Eigen<FloatType, N>&>(
+      buffer);
 }
 
-template<typename TensorType>
-auto trace(const TensorType& t) {
-    Eigen::Tensor<double, 0, Eigen::RowMajor> trace = t.value().trace();
+template<typename FloatType, unsigned short Rank>
+auto trace(const tensorwrapper::buffer::Eigen<FloatType, Rank>& t) {
+    Eigen::Tensor<FloatType, 0, Eigen::RowMajor> trace = t.value().trace();
     return trace.coeff();
 }
 
-template<typename TensorType>
-auto norm(const TensorType& t) {
-    Eigen::Tensor<double, 0, Eigen::RowMajor> norm =
+template<typename FloatType, unsigned short Rank>
+auto norm(const tensorwrapper::buffer::Eigen<FloatType, Rank>& t) {
+    Eigen::Tensor<FloatType, 0, Eigen::RowMajor> norm =
       t.value().square().sum().sqrt();
     return norm.coeff();
 }
-
 } // namespace test