Adds a Module that can normalize primitives

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_/make_libint_basis_set.hpp

@@ -26,9 +26,13 @@ namespace integrals::libint::detail_ {
 /** @brief Converts an NWX basis set object to a LibInt2 basis set object.
  *
  *  @param[in] bs The NWX basis set to be converted
+ *  @param[in] embed_normalization If true (default), normalization factors are
+ *             embedded into the contraction coefficients by libint. If false,
+ *             coefficients are used as-is without renormalization.
  *  @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 +82,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);
     }

src/integrals/libint/libint.cpp

@@ -197,6 +197,8 @@ void set_defaults(pluginplay::ModuleManager& mm) {
                      "Black Box Primitive Pair Estimator");
     mm.change_submod("CauchySchwarz Estimator", "Decontract Basis Set",
                      "Decontract Basis Set");
+    mm.change_submod("Primitive Normalization", "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");
@@ -221,6 +223,7 @@ void load_modules(pluginplay::ModuleManager& mm) {
     mm.add_module<CauchySchwarzPrimitiveEstimator>("CauchySchwarz Estimator");
     mm.add_module<PrimitiveErrorModel>("Primitive Error Model");
     mm.add_module<AnalyticError>("Analytic Error");
+    mm.add_module<PrimitiveNormalization>("Primitive Normalization");
 }
 
 #undef LOAD_LIBINT

src/integrals/libint/libint.hpp

@@ -35,6 +35,7 @@ DECLARE_MODULE(BlackBoxPrimitiveEstimator);
 DECLARE_MODULE(CauchySchwarzPrimitiveEstimator);
 DECLARE_MODULE(PrimitiveErrorModel);
 DECLARE_MODULE(AnalyticError);
+DECLARE_MODULE(PrimitiveNormalization);
 
 using simde::type::braket;
 

src/integrals/libint/primitive_normalization.cpp

@@ -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.
+ */
+
+#include "detail_/make_libint_basis_set.hpp"
+#include "libint.hpp"
+#include <cmath>
+#include <integrals/property_types.hpp>
+
+namespace integrals::libint {
+namespace {
+
+const auto desc = R"(
+Primitive Normalization
+=======================
+
+Computes the normalization factor :math:`1/\sqrt{(p|p)}` for each primitive
+in the provided basis set, where :math:`(p|p)` is the diagonal overlap
+integral of the primitive with itself.
+
+The basis set is first decontracted so that each primitive is treated as a
+separate shell. A libint2 overlap engine is then constructed with
+``embed_normalization_into_coefficients = false`` so that the raw Gaussian
+exponents and coefficients are used without any renormalization. The diagonal
+overlap :math:`(p|p)` is extracted for each primitive shell and its
+reciprocal square root is returned.
+)";
+
+} // namespace
+
+using decontract_pt = integrals::property_types::DecontractBasisSet;
+using pt = integrals::property_types::Normalize<simde::type::ao_basis_set>;
+
+MODULE_CTOR(PrimitiveNormalization) {
+    satisfies_property_type<pt>();
+    description(desc);
+
+    add_submodule<decontract_pt>("Decontract Basis Set")
+      .set_description("Module used to decontract the basis set into "
+                       "individual primitives");
+}
+
+MODULE_RUN(PrimitiveNormalization) {
+    const auto& [basis] = pt::unwrap_inputs(inputs);
+
+    // Decontract so each primitive is its own shell
+    auto& decontract_mod   = submods.at("Decontract Basis Set");
+    const auto& prim_basis = decontract_mod.run_as<decontract_pt>(basis);
+
+    // Build a libint basis set without embedding normalization into
+    // coefficients
+    auto libint_bs = detail_::make_libint_basis_set(prim_basis, false);
+
+    const auto n_shells = libint_bs.size();
+
+    // Build an overlap engine over the decontracted basis
+    if(!libint2::initialized()) libint2::initialize();
+    const auto max_nprims = libint2::max_nprim(libint_bs);
+    const auto max_l      = libint2::max_l(libint_bs);
+    libint2::Engine engine(libint2::Operator::overlap, max_nprims, max_l);
+    engine.set_max_nprim(max_nprims);
+    engine.set(libint2::BraKet::xs_xs);
+
+    const auto& buf = engine.results();
+
+    std::vector<double> norms;
+    norms.reserve(prim_basis.n_primitives());
+
+    for(std::size_t i = 0; i < n_shells; ++i) {
+        engine.compute(libint_bs[i], libint_bs[i]);
+        const auto* ints = buf[0];
+
+        // Each decontracted shell has size() AOs; extract the diagonal elements
+        const auto n_aos = libint_bs[i].size();
+        for(std::size_t mu = 0; mu < n_aos; ++mu) {
+            // Diagonal element of the (n_aos x n_aos) overlap block
+            const auto diag = ints[mu * n_aos + mu];
+            norms.push_back(1.0 / std::sqrt(diag));
+        }
+    }
+
+    auto result = results();
+    return pt::wrap_results(result, norms);
+}
+
+} // namespace integrals::libint

tests/cxx/unit/integrals/libint/detail_/test_make_libint_basis_set.cpp

@@ -24,8 +24,37 @@ using namespace integrals::testing;
 
 TEST_CASE("make_libint_basis_set") {
     using integrals::libint::detail_::make_libint_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);
+    auto aobs = water_sto3g_basis_set();
+
+    SECTION("embed_normalization=true (default)") {
+        auto libint_bs   = make_libint_basis_set(aobs);
+        auto libint_corr = test::water_basis_set();
+        REQUIRE(libint_bs == libint_corr);
+    }
+
+    SECTION("embed_normalization=false") {
+        auto libint_norm   = make_libint_basis_set(aobs, true);
+        auto libint_nonorm = make_libint_basis_set(aobs, false);
+
+        // The two basis sets should have the same number of shells
+        REQUIRE(libint_norm.size() == libint_nonorm.size());
+
+        // With normalization disabled the coefficients must differ from the
+        // renormalized ones for at least one shell (the p shell on oxygen
+        // is the clearest case since its normalization factor != 1)
+        bool any_differ = false;
+        for(std::size_t s = 0; s < libint_norm.size(); ++s) {
+            const auto& c_norm   = libint_norm[s].contr[0].coeff;
+            const auto& c_nonorm = libint_nonorm[s].contr[0].coeff;
+            REQUIRE(c_norm.size() == c_nonorm.size());
+            for(std::size_t p = 0; p < c_norm.size(); ++p) {
+                if(c_norm[p] != c_nonorm[p]) {
+                    any_differ = true;
+                    break;
+                }
+            }
+            if(any_differ) break;
+        }
+        REQUIRE(any_differ);
+    }
 }

tests/cxx/unit/integrals/libint/primitive_normalization.cpp

@@ -0,0 +1,87 @@
+/*
+ * 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 <integrals/property_types.hpp>
+
+using namespace integrals::testing;
+
+using pt = integrals::property_types::Normalize<simde::type::ao_basis_set>;
+
+namespace {
+
+/* Reference normalization factors 1/sqrt((p|p)) for water STO-3G, computed
+ * analytically from the Gaussian overlap formula.
+ *
+ * The decontracted basis order follows DecontractBasisSet:
+ *   O s1 : 3 primitives (zeta = 130.7093200, 23.8088610, 6.4436083)
+ *   O s2 : 3 primitives (zeta =   5.0331513,  1.1695961, 0.3803890)
+ *   O p  : 3 primitives × 3 AOs (zeta = 5.0331513, 1.1695961, 0.3803890)
+ *   H1 s : 3 primitives (zeta = 3.425250914, 0.6239137298, 0.1688554040)
+ *   H2 s : 3 primitives (same as H1)
+ *
+ * For s-type: (p|p) = (pi/(2*zeta))^(3/2)
+ * For p-type (e.g. x*exp(-zeta*r^2)):
+ *   (p|p) = integral x^2 exp(-2*zeta*r^2) d^3r = (1/(4*zeta)) *
+ * (pi/(2*zeta))^(3/2) Normalization factor = 1/sqrt((p|p))
+ */
+std::vector<double> corr_water_norms() {
+    return {
+      // O s1
+      27.551167600757328,
+      7.681818767185153,
+      2.882417868807197,
+      // O s2
+      2.394914875721071,
+      0.801561825779394,
+      0.345208161163625,
+      // O p (3 primitives × 3 AOs = 9 entries, same value per AO)
+      10.745832583524919,
+      10.745832583524919,
+      10.745832583524919,
+      1.733744024405248,
+      1.733744024405248,
+      1.733744024405248,
+      0.425818989418287,
+      0.425818989418287,
+      0.425818989418287,
+      // H1 s
+      1.794441833790094,
+      0.500326492211116,
+      0.187735461846361,
+      // H2 s
+      1.794441833790094,
+      0.500326492211116,
+      0.187735461846361,
+    };
+}
+
+} // namespace
+
+TEST_CASE("PrimitiveNormalization") {
+    auto mm   = initialize_integrals();
+    auto& mod = mm.at("Primitive Normalization");
+
+    SECTION("H2O/STO-3G") {
+        auto aobs  = water_sto3g_basis_set();
+        auto norms = mod.run_as<pt>(aobs);
+        auto corr  = corr_water_norms();
+
+        REQUIRE(norms.size() == corr.size());
+        for(std::size_t i = 0; i < norms.size(); ++i) {
+            REQUIRE(norms[i] == Catch::Approx(corr[i]).epsilon(1.0e-10));
+        }
+    }
+}