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_/fill_tensor.hpp b/src/integrals/libint/detail_/fill_tensor.hpp new file mode 100644 index 00000000..af596fde --- /dev/null +++ b/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 +#ifdef _OPENMP +#include +#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 +auto build_eigen_buffer(const std::vector& basis_sets, + parallelzone::runtime::RuntimeView& rv, 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); + } + auto N = basis_sets.size(); + std::vector 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(s, initial_value); + return std::make_unique(std::move(buffer)); +} + +template +auto fill_tensor(const std::vector& basis_sets, + const chemist::qm_operator::OperatorBase& op, + parallelzone::runtime::RuntimeView& rv, double thresh) { + using size_type = decltype(N); + + // Dimensional information + std::vector 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 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(basis_sets, rv, thresh); + auto data = pbuffer->get_mutable_data(); + auto span = wtf::buffer::contiguous_buffer_cast(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 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()); + + 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 +auto fill_tensor(const std::vector& 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(basis_sets, op, rv, + thresh); + } else { + throw std::runtime_error("Sigma support not enabled!"); + } + } else { + t = fill_tensor(basis_sets, op, rv, thresh); + } + return t; +} + +} // namespace integrals::libint::detail_ diff --git a/src/integrals/libint/detail_/get_basis_sets.hpp b/src/integrals/libint/detail_/get_basis_sets.hpp index d3e96f51..a08af22a 100644 --- a/src/integrals/libint/detail_/get_basis_sets.hpp +++ b/src/integrals/libint/detail_/get_basis_sets.hpp @@ -55,24 +55,31 @@ constexpr int get_n(const BraType& bra, const KetType& ket) { */ template std::vector 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 basis_sets; if constexpr(std::is_same_v) { - 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) { - 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) { - 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) { - 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; diff --git a/src/integrals/libint/detail_/make_libint_basis_set.hpp b/src/integrals/libint/detail_/make_libint_basis_set.hpp index d9b7fdf4..41bda0ed 100644 --- a/src/integrals/libint/detail_/make_libint_basis_set.hpp +++ b/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); } diff --git a/src/integrals/libint/libint.cpp b/src/integrals/libint/libint.cpp index 57fadc32..b9bede66 100644 --- a/src/integrals/libint/libint.cpp +++ b/src/integrals/libint/libint.cpp @@ -13,122 +13,11 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#include "../uncertain_types.hpp" + +#include "detail_/fill_tensor.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 "libint.hpp" -#include "libint_visitor.hpp" -#include - -#ifdef _OPENMP -#include -#endif - namespace integrals::libint { -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 -auto build_eigen_buffer(const std::vector& basis_sets, - parallelzone::runtime::RuntimeView& rv, 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); - } - auto N = basis_sets.size(); - std::vector 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(s, initial_value); - return std::make_unique(std::move(buffer)); -} - -template -auto fill_tensor(const std::vector& basis_sets, - const chemist::qm_operator::OperatorBase& op, - parallelzone::runtime::RuntimeView& rv, double thresh) { - using size_type = decltype(N); - - // Dimensional information - std::vector 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 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(basis_sets, rv, thresh); - auto data = pbuffer->get_mutable_data(); - auto span = wtf::buffer::contiguous_buffer_cast(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 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()); - - 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 TEMPLATED_MODULE_CTOR(Libint, BraKetType) { @@ -159,19 +48,7 @@ TEMPLATED_MODULE_RUN(Libint, BraKetType) { // 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, rv, - thresh); - } else { - throw std::runtime_error("Sigma support not enabled!"); - } - } else { - t = fill_tensor(basis_sets, op, rv, thresh); - } - + auto t = detail_::fill_tensor(basis_sets, op, rv, thresh, with_uq); auto result = results(); return my_pt::wrap_results(result, t); } @@ -201,6 +78,12 @@ void set_defaults(pluginplay::ModuleManager& mm) { mm.change_input("Benchmark ERI4", "Threshold", 1.0E-16); mm.change_submod("CauchySchwarz Estimator", "ERI4", "Benchmark ERI4"); mm.change_submod("Analytic Error", "ERI4s", "Benchmark ERI4"); + mm.change_submod("Raw Primitive ERI4", "Decontract Basis Set", + "Decontract Basis Set"); + mm.change_submod("Primitive Contractor ERI4", "Raw Primitive ERI4", + "Raw Primitive ERI4"); + mm.change_submod("Primitive Contractor ERI4", "Primitive Normalization", + "Primitive Normalization"); } #define LOAD_LIBINT(bra, op, ket, key) mm.add_module(key) @@ -221,6 +104,9 @@ 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("Raw Primitive ERI4"); + mm.add_module("Primitive Normalization"); + mm.add_module("Primitive Contractor ERI4"); } #undef LOAD_LIBINT diff --git a/src/integrals/libint/libint.hpp b/src/integrals/libint/libint.hpp index 46fcfa04..ef115277 100644 --- a/src/integrals/libint/libint.hpp +++ b/src/integrals/libint/libint.hpp @@ -35,6 +35,9 @@ DECLARE_MODULE(BlackBoxPrimitiveEstimator); DECLARE_MODULE(CauchySchwarzPrimitiveEstimator); DECLARE_MODULE(PrimitiveErrorModel); DECLARE_MODULE(AnalyticError); +DECLARE_MODULE(PrimitiveNormalization); +DECLARE_MODULE(RawPrimitiveERIs); +DECLARE_MODULE(PrimitiveContractor); using simde::type::braket; diff --git a/src/integrals/libint/libint_visitor.hpp b/src/integrals/libint/libint_visitor.hpp index b6d55f4d..5be3f000 100644 --- a/src/integrals/libint/libint_visitor.hpp +++ b/src/integrals/libint/libint_visitor.hpp @@ -14,8 +14,9 @@ * limitations under the License. */ #pragma once +#include "detail_/make_engine.hpp" +#include #include - namespace integrals::libint { class LibintVisitor : public chemist::qm_operator::OperatorVisitor { diff --git a/src/integrals/libint/primitive_contractor.cpp b/src/integrals/libint/primitive_contractor.cpp new file mode 100644 index 00000000..3467c105 --- /dev/null +++ b/src/integrals/libint/primitive_contractor.cpp @@ -0,0 +1,176 @@ +/* + * 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 "libint.hpp" +#include +#include + +namespace integrals::libint { +namespace { + +const auto desc = R"( +Primitive Contractor +==================== + +This module computes four-center electron repulsion integrals (ERIs) in the +contracted AO basis by explicitly contracting primitive integrals with +renormalized contraction coefficients. + +The algorithm proceeds in three steps: + +1. The "Raw Primitive ERI4" submodule is called with the original braket. + It internally decontracts each basis set (one shell per primitive, + coefficient = 1.0) and returns the raw primitive ERI tensor with libint's + automatic shell normalization disabled. The tensor has dimensions + [n_prim_aos, n_prim_aos, n_prim_aos, n_prim_aos], where n_prim_aos counts + all AO components across all decontracted shells (e.g. a p-shell with 3 + primitives contributes 9 decontracted AO indices: 3 primitives x 3 + components). + +2. The "Primitive Normalization" submodule is called on each of the four + contracted basis sets. It returns a vector of renormalized contraction + coefficients c[i], one per decontracted AO index, in the same ordering as + the decontracted basis (primitive index varying fastest within each + (shell, ao_component) block). + +3. The contracted AO integrals are formed by: + + AO[mu, nu, lam, sig] = + sum_{i in prims(mu)} sum_{j in prims(nu)} + sum_{k in prims(lam)} sum_{l in prims(sig)} + c0[i] * c1[j] * c2[k] * c3[l] * prim_ERI[i, j, k, l] + + where prims(mu) is the set of decontracted AO indices that belong to the + same contracted AO mu (same shell, same AO component, all primitives). +)"; + +// Build a map from decontracted AO index -> contracted AO index. +// For contracted shell s with n_prims primitives and n_aos AO components, +// the decontracted indices (prim p, ao_component a) all map to the contracted +// AO index first_ao[s] + a. +// NOTE: bs.size() returns the number of centers, not shells. +// bs.n_shells() returns the total number of shells (flattened). +std::vector build_prim_to_ao_map( + const simde::type::ao_basis_set& bs) { + std::vector map; + std::size_t contracted_ao = 0; + for(std::size_t s = 0; s < bs.n_shells(); ++s) { + const auto& shell = bs.shell(s); + const auto n_prims = shell.n_primitives(); + const auto n_aos = shell.size(); + for(std::size_t p = 0; p < n_prims; ++p) { + for(std::size_t a = 0; a < n_aos; ++a) { + map.push_back(contracted_ao + a); + } + } + contracted_ao += n_aos; + } + return map; +} + +} // namespace + +using my_pt = simde::ERI4; +using norm_pt = integrals::property_types::Normalize; + +MODULE_CTOR(PrimitiveContractor) { + satisfies_property_type(); + description(desc); + add_submodule("Raw Primitive ERI4"); + add_submodule("Primitive Normalization"); +} + +MODULE_RUN(PrimitiveContractor) { + const auto& [braket] = my_pt::unwrap_inputs(inputs); + + auto bra = braket.bra(); + auto ket = braket.ket(); + + const auto& bs0 = bra.first.ao_basis_set(); + const auto& bs1 = bra.second.ao_basis_set(); + const auto& bs2 = ket.first.ao_basis_set(); + const auto& bs3 = ket.second.ao_basis_set(); + + // Step 1: get raw primitive ERIs (decontracted, unnormalized). + auto& raw_mod = submods.at("Raw Primitive ERI4"); + const auto& prim_T = raw_mod.run_as(braket); + + // Step 2: get renormalized contraction coefficients for each basis set. + auto& norm_mod = submods.at("Primitive Normalization"); + const auto& c0 = norm_mod.run_as(bs0); + const auto& c1 = norm_mod.run_as(bs1); + const auto& c2 = norm_mod.run_as(bs2); + const auto& c3 = norm_mod.run_as(bs3); + + // Build primitive -> contracted AO index maps + auto map0 = build_prim_to_ao_map(bs0); + auto map1 = build_prim_to_ao_map(bs1); + auto map2 = build_prim_to_ao_map(bs2); + auto map3 = build_prim_to_ao_map(bs3); + + const Eigen::Index n0 = bs0.n_aos(); + const Eigen::Index n1 = bs1.n_aos(); + const Eigen::Index n2 = bs2.n_aos(); + const Eigen::Index n3 = bs3.n_aos(); + + const Eigen::Index np0 = c0.size(); + const Eigen::Index np1 = c1.size(); + const Eigen::Index np2 = c2.size(); + const Eigen::Index np3 = c3.size(); + + // Get raw data from the primitive tensor via Eigen TensorMap + using namespace tensorwrapper; + const auto pdata = buffer::get_raw_data(prim_T.buffer()); + + using prim_map_t = + Eigen::TensorMap>; + prim_map_t prim(pdata.data(), np0, np1, np2, np3); + + // Step 3: contract into AO basis + std::vector ao_data(n0 * n1 * n2 * n3, 0.0); + using ao_map_t = + Eigen::TensorMap>; + ao_map_t ao(ao_data.data(), n0, n1, n2, n3); + + for(Eigen::Index i = 0; i < np0; ++i) { + const double ci = c0[i]; + const Eigen::Index mu = map0[i]; + for(Eigen::Index j = 0; j < np1; ++j) { + const double cij = ci * c1[j]; + const Eigen::Index nu = map1[j]; + for(Eigen::Index k = 0; k < np2; ++k) { + const double cijk = cij * c2[k]; + const Eigen::Index lam = map2[k]; + for(Eigen::Index l = 0; l < np3; ++l) { + const Eigen::Index sig = map3[l]; + ao(mu, nu, lam, sig) += cijk * c3[l] * prim(i, j, k, l); + } + } + } + } + + // Wrap result in a TensorWrapper tensor + tensorwrapper::shape::Smooth shape( + {static_cast(n0), static_cast(n1), + static_cast(n2), static_cast(n3)}); + tensorwrapper::buffer::Contiguous buf(std::move(ao_data), shape); + simde::type::tensor t(shape, std::move(buf)); + + auto result = results(); + return my_pt::wrap_results(result, t); +} + +} // namespace integrals::libint diff --git a/src/integrals/libint/primitive_normalization.cpp b/src/integrals/libint/primitive_normalization.cpp new file mode 100644 index 00000000..e7423235 --- /dev/null +++ b/src/integrals/libint/primitive_normalization.cpp @@ -0,0 +1,92 @@ +/* + * 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 + +namespace integrals::libint { +namespace { + +const auto desc = R"( + Primitive Normalization + ======================= + + Returns the renormalized contraction coefficients for each primitive in the + provided basis set, matching libint's internal ``renorm()`` convention. + + For each contracted shell, libint applies two normalization steps: + + 1. A per-primitive factor :math:`N_p = \sqrt{2^l (2\zeta_p)^{l+3/2} / + (\sqrt{\pi^3} (2l-1)!!)}` is multiplied into each raw coefficient + :math:`d_p`. + + 2. All scaled coefficients are then divided by :math:`\sqrt{\langle\phi|\phi\rangle}`, + where :math:`\langle\phi|\phi\rangle` is the self-overlap of the contracted + shell computed with the already-scaled coefficients, so that the contracted + shell has unit norm. + + This module returns the resulting values :math:`d_p N_p / \sqrt{\langle\phi|\phi\rangle}` + by constructing the basis set via libint with ``embed_normalization = true`` + and reading the coefficients back from the resulting shell objects. + + The output vector has one entry per (primitive, AO component) pair, in the + same order as the decontracted basis: for each contracted shell, the + :math:`n_{\rm prims} \times n_{\rm AOs}` entries are listed with the primitive + index varying fastest. + )"; + +} // namespace + +using pt = integrals::property_types::Normalize; + +MODULE_CTOR(PrimitiveNormalization) { + satisfies_property_type(); + description(desc); +} + +MODULE_RUN(PrimitiveNormalization) { + const auto& [basis] = pt::unwrap_inputs(inputs); + + // Build a libint basis set with normalization embedded into coefficients. + // This triggers libint's renorm(), which applies both the per-primitive + // normalization factor and the contracted-shell unit-norm scaling. + auto libint_bs = detail_::make_libint_basis_set(basis, true); + + std::vector norms; + + for(std::size_t s = 0; s < libint_bs.size(); ++s) { + const auto& shell = libint_bs[s]; + const auto n_prims = shell.nprim(); + const auto& nwx_shell = basis.shell(s); + const auto n_aos = nwx_shell.size(); + + // For each primitive, emit one entry per AO component. All AO + // components within a primitive share the same renormalized coefficient + // (libint stores one coefficient per primitive per contraction). + for(std::size_t p = 0; p < n_prims; ++p) { + // shell.contr[0].coeff[p] holds the renormalized coefficient after + // libint's renorm(): d_p * N_p / sqrt(contracted-shell norm) + const auto c = shell.contr[0].coeff[p]; + for(std::size_t ao = 0; ao < n_aos; ++ao) { norms.push_back(c); } + } + } + + auto result = results(); + return pt::wrap_results(result, norms); +} + +} // namespace integrals::libint diff --git a/src/integrals/libint/raw_primitive_eris.cpp b/src/integrals/libint/raw_primitive_eris.cpp new file mode 100644 index 00000000..5133d3dd --- /dev/null +++ b/src/integrals/libint/raw_primitive_eris.cpp @@ -0,0 +1,105 @@ +/* + * 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_/fill_tensor.hpp" +#include "detail_/get_basis_sets.hpp" +#include "libint.hpp" +#include +#include +namespace integrals::libint { +namespace { + +const auto desc = R"( +Raw Primitive ERI4 +================== + +This module computes four-center electron repulsion integrals (ERIs) in a +fully decontracted (primitive) basis set with no libint normalization applied +to the shells. + +The input braket's bra and ket basis sets are each decontracted via the +"Decontract Basis Set" submodule, which replaces every contracted shell with +one shell per primitive (coefficient set to 1.0). The resulting decontracted +basis sets are then passed to the "ERI4" submodule with libint's automatic +shell-coefficient normalization disabled, so the raw primitive integrals are +returned without any rescaling. + +N.B. libint's normalization flag is a global static, so it is restored to +true immediately after the ERI4 call to avoid affecting other modules. +)"; + +} // namespace + +using my_pt = simde::ERI4; +using decontract_pt = integrals::property_types::DecontractBasisSet; + +MODULE_CTOR(RawPrimitiveERIs) { + satisfies_property_type(); + description(desc); + add_submodule("Decontract Basis Set"); + 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); +} + +MODULE_RUN(RawPrimitiveERIs) { + 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(); + + auto& decontract_mod = submods.at("Decontract Basis Set"); + + // Decontract each of the four basis sets + const auto& dc_bra0_bs = + decontract_mod.run_as(bra.first.ao_basis_set()); + const auto& dc_bra1_bs = + decontract_mod.run_as(bra.second.ao_basis_set()); + const auto& dc_ket0_bs = + decontract_mod.run_as(ket.first.ao_basis_set()); + const auto& dc_ket1_bs = + decontract_mod.run_as(ket.second.ao_basis_set()); + + // Wrap the decontracted ao_basis_set objects back into aos/aos_squared + simde::type::aos dc_bra0(dc_bra0_bs); + simde::type::aos dc_bra1(dc_bra1_bs); + simde::type::aos dc_ket0(dc_ket0_bs); + simde::type::aos dc_ket1(dc_ket1_bs); + + simde::type::aos_squared dc_bra(dc_bra0, dc_bra1); + simde::type::aos_squared dc_ket(dc_ket0, dc_ket1); + + // Disable libint's automatic shell normalization, compute, then restore + auto original_normalization = + libint2::Shell::do_enforce_unit_normalization(); + libint2::Shell::do_enforce_unit_normalization(false); + auto& rv = this->get_runtime(); + + auto basis_sets = detail_::get_basis_sets(dc_bra, dc_ket, false); + auto t = detail_::fill_tensor<4>(basis_sets, op, rv, thresh, with_uq); + libint2::Shell::do_enforce_unit_normalization(original_normalization); + + auto result = results(); + return my_pt::wrap_results(result, t); +} + +} // namespace integrals::libint diff --git a/tests/cxx/unit/integrals/libint/test_primitive_contractor.cpp b/tests/cxx/unit/integrals/libint/test_primitive_contractor.cpp new file mode 100644 index 00000000..a0ff7686 --- /dev/null +++ b/tests/cxx/unit/integrals/libint/test_primitive_contractor.cpp @@ -0,0 +1,50 @@ +/* + * 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" + +using namespace integrals::testing; +using tensorwrapper::operations::approximately_equal; + +TEST_CASE("Primitive Contractor ERI4") { + using test_pt = simde::ERI4; + + auto mm = initialize_integrals(); + REQUIRE(mm.count("Primitive Contractor ERI4")); + REQUIRE(mm.count("ERI4")); + + auto aobs = water_sto3g_basis_set(); + for(const auto& type : + {chemist::ShellType::cartesian, chemist::ShellType::pure}) { + for(const auto& l : {1, 2, 3}) { + aobs.shell(2).l() = l; + aobs.shell(2).pure() = type; + simde::type::aos aos(aobs); + simde::type::aos_squared aos_squared(aos, aos); + + simde::type::v_ee_type op{}; + + chemist::braket::BraKet braket(aos_squared, op, aos_squared); + auto corr_mod = mm.at("ERI4"); + auto test_mod = mm.at("Primitive Contractor ERI4"); + + auto T_contracted = test_mod.run_as(braket); + auto T_eri4 = corr_mod.run_as(braket); + + REQUIRE(approximately_equal(T_contracted, T_eri4, 1.0e-14)); + } + } +}