diff --git a/include/integrals/property_types.hpp b/include/integrals/property_types.hpp index 3998de81..f12e060e 100644 --- a/include/integrals/property_types.hpp +++ b/include/integrals/property_types.hpp @@ -91,4 +91,24 @@ TEMPLATED_PROPERTY_TYPE_RESULTS(Uncertainty, BasePT) { using DecontractBasisSet = simde::Convert; +template +DECLARE_TEMPLATED_PROPERTY_TYPE(Normalize, T); + +template +TEMPLATED_PROPERTY_TYPE_INPUTS(Normalize, T) { + auto rv = + pluginplay::declare_input().add_field("Object to Normalize"); + rv["Object to Normalize"].set_description("The object to normalize"); + return rv; +} + +template +TEMPLATED_PROPERTY_TYPE_RESULTS(Normalize, T) { + auto rv = pluginplay::declare_result().add_field>( + "Normalization Factors"); + rv["Normalization Factors"].set_description( + "A vector of normalization factors, one per primitive"); + return rv; +} + } // end namespace integrals::property_types diff --git a/src/integrals/libint/detail_/make_libint_basis_set.hpp b/src/integrals/libint/detail_/make_libint_basis_set.hpp index d9b7fdf4..21ad5ff4 100644 --- a/src/integrals/libint/detail_/make_libint_basis_set.hpp +++ b/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); } diff --git a/src/integrals/libint/libint.cpp b/src/integrals/libint/libint.cpp index 57fadc32..5c59838c 100644 --- a/src/integrals/libint/libint.cpp +++ b/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("CauchySchwarz Estimator"); mm.add_module("Primitive Error Model"); mm.add_module("Analytic Error"); + mm.add_module("Primitive Normalization"); } #undef LOAD_LIBINT diff --git a/src/integrals/libint/libint.hpp b/src/integrals/libint/libint.hpp index 46fcfa04..1455f795 100644 --- a/src/integrals/libint/libint.hpp +++ b/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; diff --git a/src/integrals/libint/primitive_normalization.cpp b/src/integrals/libint/primitive_normalization.cpp new file mode 100644 index 00000000..915e2123 --- /dev/null +++ b/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 +#include + +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; + +MODULE_CTOR(PrimitiveNormalization) { + satisfies_property_type(); + description(desc); + + add_submodule("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(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 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 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 1e4d2b57..5437cbf8 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 @@ -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); + } } diff --git a/tests/cxx/unit/integrals/libint/primitive_normalization.cpp b/tests/cxx/unit/integrals/libint/primitive_normalization.cpp new file mode 100644 index 00000000..1ebc478c --- /dev/null +++ b/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 + +using namespace integrals::testing; + +using pt = integrals::property_types::Normalize; + +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 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(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)); + } + } +}