Build ERIs from primitive quartets

include/integrals/property_types.hpp

@@ -91,4 +91,24 @@ TEMPLATED_PROPERTY_TYPE_RESULTS(Uncertainty, BasePT) {
 using DecontractBasisSet =
   simde::Convert<simde::type::ao_basis_set, simde::type::ao_basis_set>;
 
+template<typename T>
+DECLARE_TEMPLATED_PROPERTY_TYPE(Normalize, T);
+
+template<typename T>
+TEMPLATED_PROPERTY_TYPE_INPUTS(Normalize, T) {
+    auto rv =
+      pluginplay::declare_input().add_field<const T&>("Object to Normalize");
+    rv["Object to Normalize"].set_description("The object to normalize");
+    return rv;
+}
+
+template<typename T>
+TEMPLATED_PROPERTY_TYPE_RESULTS(Normalize, T) {
+    auto rv = pluginplay::declare_result().add_field<std::vector<double>>(
+      "Normalization Factors");
+    rv["Normalization Factors"].set_description(
+      "A vector of normalization factors, one per primitive");
+    return rv;
+}
+
 } // end namespace integrals::property_types

src/integrals/libint/detail_/fill_tensor.hpp

@@ -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.
+ */
+
+#pragma once
+#include "../../uncertain_types.hpp"
+#include "../libint_visitor.hpp"
+#include "libint_op.hpp"
+#include "make_engine.hpp"
+#include "make_libint_basis_set.hpp"
+#include "shells2ord.hpp"
+#include <type_traits>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+namespace integrals::libint::detail_ {
+namespace {
+#ifdef _OPENMP
+int get_num_threads() {
+    int num_threads;
+#pragma omp parallel
+    { num_threads = omp_get_num_threads(); }
+    return num_threads;
+}
+
+int get_thread_num() { return omp_get_thread_num(); }
+#else
+
+int get_num_threads() { return 1; }
+
+int get_thread_num() { return 0; }
+
+#endif
+
+template<typename FloatType>
+auto build_eigen_buffer(const std::vector<libint2::BasisSet>& basis_sets,
+                        parallelzone::runtime::RuntimeView& rv, 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);
+    }
+    auto N = basis_sets.size();
+    std::vector<decltype(N)> dims(N);
+    for(decltype(N) i = 0; i < N; ++i) dims[i] = basis_sets[i].nbf();
+
+    using namespace tensorwrapper;
+    using shape_t = shape::Smooth;
+
+    shape_t s{dims.begin(), dims.end()};
+    auto buffer = buffer::make_contiguous<FloatType>(s, initial_value);
+    return std::make_unique<buffer::Contiguous>(std::move(buffer));
+}
+
+template<std::size_t N, typename FloatType>
+auto fill_tensor(const std::vector<libint2::BasisSet>& basis_sets,
+                 const chemist::qm_operator::OperatorBase& op,
+                 parallelzone::runtime::RuntimeView& rv, double thresh) {
+    using size_type = decltype(N);
+
+    // Dimensional information
+    std::vector<size_type> dim_stepsizes(N, 1);
+    size_type num_shell_combinations = 1;
+
+    for(size_type i = 0; i < N; ++i) {
+        num_shell_combinations *= basis_sets[i].size();
+        for(size_type j = i; j < N - 1; ++j) {
+            dim_stepsizes[i] *= basis_sets[j].size();
+        }
+    }
+
+    // Make an engine for each thread
+    int num_threads = get_num_threads();
+    std::vector<libint2::Engine> engines(num_threads);
+    LibintVisitor visitor(basis_sets, thresh);
+    op.visit(visitor);
+    for(int i = 0; i != num_threads; ++i) { engines[i] = visitor.engine(); }
+
+    // Fill in values
+    auto pbuffer = build_eigen_buffer<FloatType>(basis_sets, rv, thresh);
+    auto data    = pbuffer->get_mutable_data();
+    auto span    = wtf::buffer::contiguous_buffer_cast<FloatType>(data);
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+    for(size_type i_pair = 0; i_pair != num_shell_combinations; ++i_pair) {
+        auto thread_id = get_thread_num();
+
+        std::vector<size_type> shells(N);
+        auto shell_ord = i_pair;
+        for(size_type i = 0; i < N; ++i) {
+            shells[i] = shell_ord / dim_stepsizes[i];
+            shell_ord = shell_ord % dim_stepsizes[i];
+        }
+
+        detail_::run_engine_(engines[thread_id], basis_sets, shells,
+                             std::make_index_sequence<N>());
+
+        const auto& buf = engines[thread_id].results();
+        auto vals       = buf[0];
+        if(vals) {
+            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) {
+                auto update      = span[ord[i_ord]] + vals[i_ord];
+                span[ord[i_ord]] = update;
+            }
+        }
+    }
+
+    auto pshape = pbuffer->layout().shape().clone();
+    return simde::type::tensor(std::move(pshape), std::move(pbuffer));
+}
+} // namespace
+
+template<std::size_t N>
+auto fill_tensor(const std::vector<libint2::BasisSet>& basis_sets,
+                 const chemist::qm_operator::OperatorBase& op,
+                 parallelzone::runtime::RuntimeView& rv, double thresh,
+                 bool with_uq) {
+    simde::type::tensor t;
+    if(with_uq) {
+        if constexpr(integrals::type::has_sigma()) {
+            t = fill_tensor<N, type::uncertain_double>(basis_sets, op, rv,
+                                                       thresh);
+        } else {
+            throw std::runtime_error("Sigma support not enabled!");
+        }
+    } else {
+        t = fill_tensor<N, double>(basis_sets, op, rv, thresh);
+    }
+    return t;
+}
+
+} // namespace integrals::libint::detail_

src/integrals/libint/detail_/get_basis_sets.hpp

@@ -55,24 +55,31 @@ constexpr int get_n(const BraType& bra, const KetType& ket) {
  */
 template<typename BraType, typename KetType>
 std::vector<libint2::BasisSet> get_basis_sets(const BraType& bra,
-                                              const KetType& ket) {
+                                              const KetType& ket,
+                                              bool embed_normalization = true) {
     using simde::type::aos;
     using simde::type::aos_squared;
 
     std::vector<libint2::BasisSet> basis_sets;
 
     if constexpr(std::is_same_v<BraType, aos>) {
-        basis_sets.push_back(make_libint_basis_set(bra.ao_basis_set()));
+        basis_sets.push_back(
+          make_libint_basis_set(bra.ao_basis_set(), embed_normalization));
     } else if constexpr(std::is_same_v<BraType, aos_squared>) {
-        basis_sets.push_back(make_libint_basis_set(bra.first.ao_basis_set()));
-        basis_sets.push_back(make_libint_basis_set(bra.second.ao_basis_set()));
+        basis_sets.push_back(
+          make_libint_basis_set(bra.first.ao_basis_set(), embed_normalization));
+        basis_sets.push_back(make_libint_basis_set(bra.second.ao_basis_set(),
+                                                   embed_normalization));
     }
 
     if constexpr(std::is_same_v<KetType, aos>) {
-        basis_sets.push_back(make_libint_basis_set(ket.ao_basis_set()));
+        basis_sets.push_back(
+          make_libint_basis_set(ket.ao_basis_set(), embed_normalization));
     } else if constexpr(std::is_same_v<KetType, aos_squared>) {
-        basis_sets.push_back(make_libint_basis_set(ket.first.ao_basis_set()));
-        basis_sets.push_back(make_libint_basis_set(ket.second.ao_basis_set()));
+        basis_sets.push_back(
+          make_libint_basis_set(ket.first.ao_basis_set(), embed_normalization));
+        basis_sets.push_back(make_libint_basis_set(ket.second.ao_basis_set(),
+                                                   embed_normalization));
     }
 
     return basis_sets;

src/integrals/libint/detail_/make_libint_basis_set.hpp

@@ -28,7 +28,8 @@ namespace integrals::libint::detail_ {
  *  @param[in] bs The NWX basis set to be converted
  *  @returns The basis set as a LibInt2 basis set
  */
-inline auto make_libint_basis_set(const simde::type::ao_basis_set& bs) {
+inline auto make_libint_basis_set(const simde::type::ao_basis_set& bs,
+                                  bool embed_normalization = true) {
     /// Typedefs for everything
     using atom_t          = libint2::Atom;
     using shell_t         = libint2::Shell;
@@ -78,7 +79,8 @@ inline auto make_libint_basis_set(const simde::type::ao_basis_set& bs) {
             svec_d_t coefs(&prim0.coefficient(), &primN.coefficient() + 1);
             conts_t conts{cont_t{l, pure, coefs}};
             /// Use origin for position, because BasisSet moves shells to center
-            atom_bases.push_back(shell_t(alphas, conts, origin));
+            atom_bases.push_back(
+              shell_t(alphas, conts, origin, embed_normalization));
         }
         element_bases.push_back(atom_bases);
     }