From 505590c0b5f00304ca3afbe3188585c1f790e9f0 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Sun, 25 Oct 2026 23:02:10 -0400 Subject: [PATCH 01/11] Engine: check the angular momentum against max_l in all builds; consistent lmax_exceeded MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - compute1()/compute2() only asserted that the shells do not exceed the max_l of the Engine, hence a Release build overran the stack and scratch buffers sized for max_l; throw Engine::lmax_exceeded instead. - lmax_exceeded thrown by initialize() reported the max AM + 1 as the limit, the one thrown by compute2() the max AM; report the max AM in both, and include the task and the limits in what(). - rkb_operators: the "folded" paragraph applies to the operators with one folded electron; opop_coulomb_opop already contains the -1 of i², as its prefactor column says. - tests for the per-center checks of compute2(), the max_l check, the limit reported by initialize(), and the core evaluator rebuilt by set(BraKet). --- export/tests/unit/test-2body.cc | 77 +++++++++++++++++++++++++++++++++ include/libint2/engine.h | 20 +++++++-- include/libint2/engine.impl.h | 12 ++--- 3 files changed, 101 insertions(+), 8 deletions(-) diff --git a/export/tests/unit/test-2body.cc b/export/tests/unit/test-2body.cc index d346135e1..621648339 100644 --- a/export/tests/unit/test-2body.cc +++ b/export/tests/unit/test-2body.cc @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -741,6 +742,82 @@ TEST_CASE("rejected set() leaves no kernels", "[engine][2-body]") { #endif } +// compute() must throw Engine::lmax_exceeded, in all builds, for a shell +// beyond the max_l of the Engine or beyond the limit of the kernels for its +// center +TEST_CASE("compute() checks the angular momentum limits", "[engine][2-body]") { + const auto eps = std::numeric_limits::epsilon(); + const Shell s{{1.0}, {{0, false, {1.0}}}, {{0.0, 0.0, 0.0}}}; + // @return the lmax_exceeded thrown by @p f + auto thrown_by = + [](const std::function &f) -> std::pair { + try { + f(); + } catch (Engine::lmax_exceeded &e) { + return {e.lmax_limit(), e.lmax_requested()}; + } + FAIL("lmax_exceeded was not thrown"); + return {0, 0}; + }; + +#if defined(LIBINT2_SUPPORT_ERI) && LIBINT2_MAX_AM_eri > 0 + SECTION("max_l of the Engine") { + const Shell p{{1.0}, {{1, false, {1.0}}}, {{0.0, 0.0, 1.0}}}; + Engine engine(Operator::coulomb, 1, 0); + CHECK(thrown_by([&]() { engine.compute(s, s, s, p); }) == + std::make_pair(size_t(0), size_t(1))); + Engine engine1(Operator::overlap, 1, 0); + CHECK(thrown_by([&]() { engine1.compute(p, s); }) == + std::make_pair(size_t(0), size_t(1))); + } +#endif + + // the kernels of a center-dependent task support higher AM on the fitting + // center than on the AO centers, the Engine's max_l does not tell them apart +#if defined(LIBINT2_SUPPORT_ERI3) && LIBINT2_CENTER_DEPENDENT_MAX_AM_3eri && \ + LIBINT2_MAX_AM_3eri > LIBINT2_MAX_AM_default + SECTION("AM limit of the AO centers of 3-center kernels") { + const int l_ao = LIBINT2_MAX_AM_default + 1; + const Shell ao{{1.0}, {{l_ao, true, {1.0}}}, {{0.0, 0.0, 1.0}}}; + Engine engine(Operator::coulomb, 1, LIBINT2_MAX_AM_3eri, 0, eps, + libint2::default_params(Operator::coulomb), BraKet::xs_xx); + CHECK(thrown_by([&]() { engine.compute(s, ao, s); }) == + std::make_pair(size_t(LIBINT2_MAX_AM_default), size_t(l_ao))); + } +#endif + + // the limit reported by initialize() is the max AM, not the max AM + 1 +#if defined(LIBINT2_SUPPORT_ERI) + SECTION("AM limit reported by initialize()") { + CHECK(thrown_by([&]() { + Engine engine(Operator::coulomb, 1, LIBINT2_MAX_AM_eri + 1); + }) == std::make_pair(size_t(LIBINT2_MAX_AM_eri), + size_t(LIBINT2_MAX_AM_eri + 1))); + } +#endif +} + +// set(BraKet) must resize the core evaluator for the new braket: a 3-center +// Engine needs Boys functions up to 3·max_l, a 4-center one up to 4·max_l +TEST_CASE("set(BraKet) rebuilds the core evaluator", "[engine][2-body]") { +#if defined(LIBINT2_SUPPORT_ERI) && defined(LIBINT2_SUPPORT_ERI3) + const int l = std::min(LIBINT2_MAX_AM_eri, LIBINT2_MAX_AM_3eri); + const Shell a{{1.0, 0.3}, {{l, false, {0.4, 0.7}}}, {{0.0, 0.0, 0.0}}}; + const Shell b{{1.1, 0.2}, {{l, false, {0.5, 0.6}}}, {{0.0, 0.5, 0.1}}}; + const auto eps = std::numeric_limits::epsilon(); + Engine engine(Operator::coulomb, 2, l, 0, eps, + libint2::default_params(Operator::coulomb), BraKet::xs_xx); + engine.set(BraKet::xx_xx); + Engine ref(Operator::coulomb, 2, l); + const auto &result = engine.compute(a, b, b, a); + const auto &result_ref = ref.compute(a, b, b, a); + REQUIRE(result[0] != nullptr); + const auto n = a.size() * b.size() * b.size() * a.size(); + for (size_t i = 0; i != n; ++i) + REQUIRE(result[0][i] == Approx(result_ref[0][i]).margin(1e-14)); +#endif +} + TEST_CASE("RKB Coulomb integrals", "[engine][2-body]") { const std::vector all_shells{ // pseudorandom s diff --git a/include/libint2/engine.h b/include/libint2/engine.h index badd60654..32e1b12b6 100644 --- a/include/libint2/engine.h +++ b/include/libint2/engine.h @@ -54,6 +54,7 @@ #include #include #include +#include #include #include #include @@ -308,8 +309,13 @@ enum class Operator { σ_i σ_j = δ_ij + i ε_ijk σ_k into a scalar (component 0) and a vector (components 1-3) part; the user assembles the 2×2 spin block of that electron as q₀·1 + i (q₁ σ_x + q₂ σ_y + q₃ σ_z) and must not apply the Pauli matrices - again. "Not folded" means the result is the spatial tensor and the user - contracts it with the Pauli matrices. + again. This applies as stated to the operators with one folded electron + (`coulomb_opop` and the 1-body operators); for `opop_coulomb_opop`, which + folds both electrons, use its prefactor column instead: the product of the two + i's of the vector ⊗ vector block (i² = −1) is already in the kernel, so + assembling each electron's block as q₀ + i σ·q would apply it twice. + "Not folded" means the result is the spatial tensor and the user contracts it + with the Pauli matrices. | Engine (C++ API) | Compiler | task label | #comp | σ·p acts on | component index | prefactor | |---------------------|-------------------|---------------------|-------|-------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------| @@ -1273,10 +1279,18 @@ class Engine { /// Exception class to be used when the angular momentum limit is exceeded. class lmax_exceeded : virtual public std::logic_error { public: + /// @param task_name the label of the task (or the empty string) + /// @param lmax_limit the max angular momentum supported + /// @param lmax_requested the angular momentum that was requested lmax_exceeded(const char* task_name, size_t lmax_limit, size_t lmax_requested) : std::logic_error( - "Engine::lmax_exceeded -- angular momentum limit exceeded"), + std::string( + "Engine::lmax_exceeded -- angular momentum limit exceeded") + + (task_name[0] != '\0' ? std::string(" for task ") + task_name + : std::string()) + + ": the max angular momentum is " + std::to_string(lmax_limit) + + ", requested " + std::to_string(lmax_requested)), lmax_limit_(lmax_limit), lmax_requested_(lmax_requested) { strncpy(task_name_, task_name, 64); diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index e1a90f254..2cd0e4a11 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -217,8 +217,9 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute1( const auto l1 = s1.contr[0].l; const auto l2 = s2.contr[0].l; - assert(l1 <= lmax_ && "the angular momentum limit is exceeded"); - assert(l2 <= lmax_ && "the angular momentum limit is exceeded"); + // the scratch and the stack are sized for lmax_, check in all builds + if (std::max(l1, l2) > lmax_) + throw lmax_exceeded(task_name_, lmax_, std::max(l1, l2)); // if want nuclear, make sure there is at least one nucleus .. otherwise the // user likely forgot to call set_params @@ -274,7 +275,6 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute1( // adjust max angular momentum, if needed const auto lmax = std::max(l1, l2); - assert(lmax <= lmax_ && "the angular momentum limit is exceeded"); // N.B. for l=0 no need to transform to solid harmonics // this is a workaround for the corner case of oper_ == Operator::*nuclear, @@ -660,7 +660,7 @@ __libint2_engine_inline void Engine::_initialize() { std::max(hard_lmax_, hard_default_lmax_), hard_lmax_); \ if (lmax_ >= lmax) { \ buildfnptrs_ = nullptr; /* leave no table for kernels we cannot use */ \ - throw Engine::lmax_exceeded(BOOST_PP_STRINGIZE(BOOST_PP_NBODYENGINE_MCR3_TASK(product)), lmax, lmax_); \ + throw Engine::lmax_exceeded(BOOST_PP_STRINGIZE(BOOST_PP_NBODYENGINE_MCR3_TASK(product)), lmax - 1, lmax_); \ } \ if (stack_size_ > 0) \ LIBINT2_PREFIXED_NAME(libint2_cleanup_default)(&primdata_[0]); \ @@ -1437,7 +1437,9 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( // adjust max angular momentum, if needed auto lmax = std::max(std::max(bra1.contr[0].l, bra2.contr[0].l), std::max(ket1.contr[0].l, ket2.contr[0].l)); - assert(lmax <= lmax_ && "the angular momentum limit is exceeded"); + // the scratch, the stack and the core evaluator are sized for lmax_, check + // in all builds + if (lmax > lmax_) throw lmax_exceeded(task_name_, lmax_, lmax); const auto lmax_bra = std::max(bra1.contr[0].l, bra2.contr[0].l); const auto lmax_ket = std::max(ket1.contr[0].l, ket2.contr[0].l); From ce6580f69d5d42b6543c3a15941e60447a147761 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Sun, 25 Oct 2026 23:02:29 -0400 Subject: [PATCH 02/11] hartree-fock++: accept a skipped derivative only if the kernels cannot do it; test the hessian in CI again - the validator accepted "2-body forces/hessian skipped" whenever it saw it; it now reads the max AM of the ERI kernels for each derivative order from the features file (new LIBINT_ERI_MAX_AM_DERIVS) and accepts a skip only if that is below the AM of the basis. - the CI lane with LIBINT2_MAX_AM="2;2;1" generated 2nd-derivative ERIs up to p only, so hartree-fock++ skipped its hessian, which then was not tested in any lane; generate them up to d (LIBINT2_ERI_MAX_AM="2;2;2"). - api_basic_compile_test: the 2-center derivative checks also need the 2-center (2eri) kernels. --- .github/workflows/cmake.yml | 9 ++++--- cmake/modules/int_am.cmake | 9 +++++++ cmake/modules/int_computed.cmake.in | 1 + .../hartree-fock/hartree-fock++-validate.py | 27 +++++++++++++------ export/tests/hartree-fock/hartree-fock++.cc | 17 +++++++++--- src/lib/libint/features.cmake.in | 1 + 6 files changed, 50 insertions(+), 14 deletions(-) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index e10e10a3c..5cf60d69a 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -45,15 +45,18 @@ jobs: build_shared_libs: ON os: macos-15 do_fortran: ON - # exercise a center-dependent AM table for the 3-center RKB ERIs and - # 2nd geometric derivatives of the 4-center ERIs in one lane only + # exercise a center-dependent AM table for the 3-center RKB ERIs, + # per-derivative AM lists, and 2nd geometric derivatives of the + # 4-center ERIs in one lane only; the 2nd-derivative ERIs must cover + # the d shells of the basis of hartree-fock++, else its hessian is + # skipped - build_type: Release os: macos-15 subproject: OFF rkb_eri3_max_am: 3 eri_deriv_order: 2 max_am: '2;2;1' - eri_max_am: '2;2;1' + eri_max_am: '2;2;2' name: "Repo • ${{ matrix.os }}: ${{ matrix.cxx }} ${{ matrix.build_type }} shared=${{ matrix.build_shared_libs }} fetch=${{ matrix.subproject }}" runs-on: ${{ matrix.os }} diff --git a/cmake/modules/int_am.cmake b/cmake/modules/int_am.cmake index cabcc18bd..4f0ee80e9 100644 --- a/cmake/modules/int_am.cmake +++ b/cmake/modules/int_am.cmake @@ -357,6 +357,15 @@ endmacro() process_integrals_class(ONEBODY) process_integrals_class(ERI) +# max AM of the 4-center ERIs for each derivative order, comma-separated +# (for the features file read by the Hartree-Fock tests) +set(Libint2_MAX_AM_ERI_DERIVS "") +if (LIBINT_INCLUDE_ERI GREATER -1) + foreach(_d RANGE ${LIBINT_INCLUDE_ERI}) + list(APPEND Libint2_MAX_AM_ERI_DERIVS ${_candidate_ERI_d${_d}}) + endforeach() +endif() +list(JOIN Libint2_MAX_AM_ERI_DERIVS "," Libint2_MAX_AM_ERI_DERIVS) process_integrals_class(RKB_ERI) process_integrals_class(RKB_ERI3) process_integrals_class(ERI3) diff --git a/cmake/modules/int_computed.cmake.in b/cmake/modules/int_computed.cmake.in index 51d79d527..4a66b89e2 100644 --- a/cmake/modules/int_computed.cmake.in +++ b/cmake/modules/int_computed.cmake.in @@ -87,3 +87,4 @@ list(PREPEND Libint2_CONFIG_COMPONENTS "${_o1st}${_o2nd}") # for features and libint2-config.cmake files set(Libint2_MAX_AM_ERI "@_candidate_ERI_d0@") +set(Libint2_MAX_AM_ERI_DERIVS "@Libint2_MAX_AM_ERI_DERIVS@") diff --git a/export/tests/hartree-fock/hartree-fock++-validate.py b/export/tests/hartree-fock/hartree-fock++-validate.py index 95be09631..8e6f2a929 100644 --- a/export/tests/hartree-fock/hartree-fock++-validate.py +++ b/export/tests/hartree-fock/hartree-fock++-validate.py @@ -44,6 +44,21 @@ def validate(label, data, refdata, tolerance, textline): if LIBINT_ERI_MAX_AM<2: sys.exit(0) +# the max AM of the basis (aug-cc-pVDZ for H2O) +OBS_MAX_L = 2 + +def skip_ok(label, deriv_order, textline): + """a quantity that needs the ERI derivatives of order deriv_order may only + be skipped if their kernels do not support the AM of the basis""" + expected = LIBINT_ERI_MAX_AM_DERIVS[deriv_order] < OBS_MAX_L + if expected: + print(textline, end="") + else: + print(label, "check: failed\nskipped although the library supports l <=", + LIBINT_ERI_MAX_AM_DERIVS[deriv_order], "for derivative order", + deriv_order, "\nactual:", textline) + return expected + eref = [-76.003354058439] etol = 5e-12 @@ -175,17 +190,13 @@ def validate(label, data, refdata, tolerance, textline): elif match15: sQok = validate("spherical electric quadrupole moment", match15.groups(), sQref, sQtol, line) elif match16: - F2ok = True - print(line,end="") + F2ok = skip_ok("2-body force", 1, line) elif match17: - Fok = True - print(line,end="") + Fok = skip_ok("HF force", 1, line) elif match18: - H2ok = True - print(line,end="") + H2ok = skip_ok("2-body hessian", 2, line) elif match19: - Hok = True - print(line,end="") + Hok = skip_ok("HF hessian", 2, line) else: print(line,end="") diff --git a/export/tests/hartree-fock/hartree-fock++.cc b/export/tests/hartree-fock/hartree-fock++.cc index 2e45d3b41..ac45d1d9c 100644 --- a/export/tests/hartree-fock/hartree-fock++.cc +++ b/export/tests/hartree-fock/hartree-fock++.cc @@ -2589,8 +2589,14 @@ void api_basic_compile_test(const BasisSet& obs, } #if LIBINT2_DERIV_ERI_ORDER - // skipped if the derivative kernels do not support the AM of the basis - if (obs.max_l() <= LIBINT2_MAX_AM_eri1) { // test deriv 2-index ints + // skipped if the derivative kernels do not support the AM of the basis; this + // needs the 4-center (eri1) and the 2-center (2eri1) kernels +#if defined(LIBINT2_MAX_AM_2eri1) + const int lmax_d1 = std::min(LIBINT2_MAX_AM_eri1, LIBINT2_MAX_AM_2eri1); +#else + const int lmax_d1 = LIBINT2_MAX_AM_eri1; +#endif + if (obs.max_l() <= lmax_d1) { // test deriv 2-index ints Engine eri4_engine(Operator::coulomb, obs.max_nprim(), obs.max_l(), 1); Engine eri2_engine = eri4_engine; eri2_engine.set(BraKet::xs_xs); @@ -2625,7 +2631,12 @@ void api_basic_compile_test(const BasisSet& obs, #endif #if LIBINT2_DERIV_ERI_ORDER > 1 - if (obs.max_l() <= LIBINT2_MAX_AM_eri2) { // test 2nd deriv 2-index ints +#if defined(LIBINT2_MAX_AM_2eri2) + const int lmax_d2 = std::min(LIBINT2_MAX_AM_eri2, LIBINT2_MAX_AM_2eri2); +#else + const int lmax_d2 = LIBINT2_MAX_AM_eri2; +#endif + if (obs.max_l() <= lmax_d2) { // test 2nd deriv 2-index ints Engine eri4_engine(Operator::coulomb, obs.max_nprim(), obs.max_l(), 2); Engine eri2_engine = eri4_engine; eri2_engine.set(BraKet::xs_xs); diff --git a/src/lib/libint/features.cmake.in b/src/lib/libint/features.cmake.in index 2d310e873..0d3fff762 100644 --- a/src/lib/libint/features.cmake.in +++ b/src/lib/libint/features.cmake.in @@ -3,6 +3,7 @@ LIBINT_ONEBODY_DERIV = @LIBINT_ONEBODY_DERIV@ LIBINT_SUPPORTS_ERI = @LIBINT_SUPPORTS_ERI@ LIBINT_ERI_MAX_AM = @Libint2_MAX_AM_ERI@ LIBINT_ERI_DERIV = @LIBINT_ERI_DERIV@ +LIBINT_ERI_MAX_AM_DERIVS = [@Libint2_MAX_AM_ERI_DERIVS@] LIBINT_CONTRACTED_INTS = @LIBINT_CONTRACTED_INTS@ LIBINT_SHELL_SET = @LIBINT_SHELL_SET@ LIBINT_SHGSHELL_ORDERING = @LIBINT_SHGSHELL_ORDERING@ From 1b09a6d028da85a965ab6417073293d2379a25bd Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Sun, 25 Oct 2026 23:02:35 -0400 Subject: [PATCH 03/11] Compiler: RKB/DerivGaussV2 follow-ups MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - rkb_fold_make_children: lead the symbolic sum with a +1 term, so that the cross x cross components of opop_coulomb_opop no longer multiply by -1.0 (which was also missing from the flop count). - descriptors of the RKB and σ·p 1-body operators: check the component index (both bounds) in all builds. - build_libint: PairSwap::coupled requires p1p2_swappable, check at compile time instead of asserting in the predicates; stale comments (translational invariance is handled by DerivGaussV2), unused includes, a leftover flush. - master_ints_list: SMultipole_1_1_int was outside the LIBINT_SUPPORT_ONEBODYINTS guard of the list. - comments: the reason for the DerivGaussV2TI label, "quaternions". --- src/bin/libint/build_libint.cc | 15 +++-- ...3p\317\203pCoulomb\317\203p\317\203p_11.h" | 8 +-- src/bin/libint/comp_deriv_gauss_v2.h | 4 +- src/bin/libint/master_ints_list.h | 9 +-- src/bin/libint/oper.h | 23 +++++--- src/bin/libint/rkb_fold_codegen.h | 58 ++++++++++++------- 6 files changed, 69 insertions(+), 48 deletions(-) diff --git a/src/bin/libint/build_libint.cc b/src/bin/libint/build_libint.cc index 2726ea585..3f95d4c5e 100644 --- a/src/bin/libint/build_libint.cc +++ b/src/bin/libint/build_libint.cc @@ -31,9 +31,7 @@ */ #include -#include #include -#include #include #include #include @@ -95,7 +93,6 @@ struct ShellQuartetSetPredicate { case PairSwap::coupled: // the ket is sorted and the bra follows; when the totals tie the // engine swaps bra and ket such that max(la,lb) <= lc - assert(p1p2_swappable); return lc >= ld && (la + lb < lc + ld || (la + lb == lc + ld && std::max(la, lb) <= lc)); case PairSwap::none: @@ -113,7 +110,6 @@ struct ShellQuartetSetPredicate { return la <= lb && lc <= ld && (!p1p2_swappable || (la < lc || (la == lc && lb <= ld))); case PairSwap::coupled: - assert(p1p2_swappable); return lc <= ld && (la + lb > lc + ld || (la + lb == lc + ld && std::min(la, lb) >= lc)); case PairSwap::none: @@ -1312,6 +1308,10 @@ static void build_TwoPRep_2b_2k( // (see the descriptor's p1p2_swappable and pair_swap) constexpr bool p1p2_swappable = OperDescrType::p1p2_swappable; constexpr PairSwap pair_swap = OperDescrType::pair_swap; + // the shell-set predicates assume that coupled pair swaps come with a + // particle swap (the engine swaps the pairs together with bra and ket) + static_assert(pair_swap != PairSwap::coupled || p1p2_swappable, + "PairSwap::coupled requires a p1p2-swappable operator"); // Note: la, lb, lc, ld generate code for chemist notation (ab|O|cd), where O // is a two-body operator. @@ -1377,7 +1377,7 @@ static void build_TwoPRep_2b_2k( //////////// // loop over unique derivative index combinations //////////// - // NB translational invariance is now handled by CR_DerivGauss + // NB translational invariance is handled by DerivGaussV2 CartesianDerivIterator<4> diter(deriv_level); bool last_deriv = false; @@ -1609,7 +1609,7 @@ static void build_TwoPRep_1b_2k( //////////// // loop over unique derivative index combinations //////////// - // NB translational invariance is now handled by CR_DerivGauss + // NB translational invariance is handled by DerivGaussV2 CartesianDerivIterator<3> diter(deriv_level); std::vector> targets; bool last_deriv = false; @@ -1832,7 +1832,7 @@ void build_TwoPRep_1b_1k(std::ostream& os, //////////// // loop over unique derivative index combinations //////////// - // NB translational invariance is now handled by CR_DerivGauss + // NB translational invariance is handled by DerivGaussV2 CartesianDerivIterator<2> diter(deriv_level); std::vector> targets; bool last_deriv = false; @@ -1909,7 +1909,6 @@ void build_TwoPRep_1b_1k(std::ostream& os, g_progress.current_task = label; g_progress.print(); - std::cout.flush(); std::string prefix(cparams->source_directory()); std::deque decl_filenames; diff --git "a/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h" "b/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h" index 0404da395..212d2d164 100644 --- "a/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h" +++ "b/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h" @@ -25,11 +25,10 @@ #include #include #include +#include #include #include -#include - #include namespace libint2 { @@ -48,8 +47,9 @@ namespace libint2 { * Molecules", Springer, 2007): * \f$ (\sigma \cdot u)(\sigma \cdot v) = (u \cdot v)I * + i\sigma \cdot (u \times v) \f$ - * applied independently to each electron's spin space gives a tensor product - * of two quaternions with \f$ 4 \times 4 = 16 \f$ components: + * applied independently to each electron's spin space gives the tensor + * product of two 4-component (scalar + vector) objects, with + * \f$ 4 \times 4 = 16 \f$ components: * * index = 4 * q₁ + q₂ * diff --git a/src/bin/libint/comp_deriv_gauss_v2.h b/src/bin/libint/comp_deriv_gauss_v2.h index 172540677..018158826 100644 --- a/src/bin/libint/comp_deriv_gauss_v2.h +++ b/src/bin/libint/comp_deriv_gauss_v2.h @@ -436,8 +436,8 @@ std::string DerivGaussV2 aux0(new mType(0u)); - // N.B. not "CR_DerivGauss...": CR_DerivGauss labels its own instances - // that way, and the RRStack deduplicates by label + // N.B. a prefix of its own, distinct from the "CR_DerivGauss " labels of + // CR_DerivGauss instances (the RRStack deduplicates by label) os << "DerivGaussV2TI" << "P" << part << to_string(where) << "D" << dir_ << genintegralset_label(target_->bra(), target_->ket(), aux0, diff --git a/src/bin/libint/master_ints_list.h b/src/bin/libint/master_ints_list.h index 8ed891e8f..af528d1d3 100644 --- a/src/bin/libint/master_ints_list.h +++ b/src/bin/libint/master_ints_list.h @@ -159,16 +159,17 @@ typedef boost::mpl::list< σpeμσp_1_1_int, CMultipole_1_1_sh, CMultipole_1_1_int, CMultipole_1_1_sh_x, CMultipole_1_1_sh_y, CMultipole_1_1_sh_z, CMultipole_1_1_int_x, CMultipole_1_1_int_y, CMultipole_1_1_int_z, SMultipole_1_1_sh, + SMultipole_1_1_int #endif - SMultipole_1_1_int> + > MasterIntegralTypeList_1body_part; typedef boost::mpl::list< TwoPRep_11_11_sq, TwoPRep_11_11_int, Coulombσpσp_11_11_sq, Coulombσpσp_11_11_int, σpσpCoulombσpσp_11_11_sq, σpσpCoulombσpσp_11_11_int, σpCoulombσp_11_11_sq, σpCoulombσp_11_11_int, Coulombσp_11_11_sq, - Coulombσp_11_11_int, R12kG12_11_11_sq, - R12kG12_11_11_int, R12kR12lG12_11_11_sq, R12kR12lG12_11_11_int, - TiG12_11_11_sq, TiG12_11_11_int, G12TiG12_11_11_sq, G12TiG12_11_11_int, + Coulombσp_11_11_int, R12kG12_11_11_sq, R12kG12_11_11_int, + R12kR12lG12_11_11_sq, R12kR12lG12_11_11_int, TiG12_11_11_sq, + TiG12_11_11_int, G12TiG12_11_11_sq, G12TiG12_11_11_int, DivG12prime_xTx_11_11_sq, DivG12prime_xTx_11_11_int, DummySymmIntegral_11_11_sq, DummySymmIntegral_11_11_int> MasterIntegralTypeList_2body_part; diff --git a/src/bin/libint/oper.h b/src/bin/libint/oper.h index aec9a8867..100f5173f 100644 --- a/src/bin/libint/oper.h +++ b/src/bin/libint/oper.h @@ -30,6 +30,7 @@ #include #include +#include #include namespace libint2 { @@ -304,7 +305,8 @@ struct σpVσp_Descr : public Contractable<σpVσp_Descr> { σpVσp_Descr() : quaternion_index_(0) {} σpVσp_Descr(int quaternion_index) : quaternion_index_(quaternion_index) { - assert(quaternion_index <= 3); + if (quaternion_index < 0 || quaternion_index > 3) + throw std::invalid_argument("quaternion_index out of range"); } static const unsigned int max_key = 4; @@ -345,7 +347,8 @@ struct σpeμσp_Descr : public Contractable<σpeμσp_Descr> { σpeμσp_Descr() : composite_index_(0) {} σpeμσp_Descr(int composite_index) : composite_index_(composite_index) { - assert(composite_index >= 0 && composite_index < 12); + if (composite_index < 0 || composite_index > 11) + throw std::invalid_argument("composite_index out of range"); } static const unsigned int max_key = 12; @@ -472,7 +475,8 @@ struct Coulombσpσp_Descr : public Contractable { Coulombσpσp_Descr() : quaternion_index_(0) {} Coulombσpσp_Descr(int quaternion_index) : quaternion_index_(quaternion_index) { - assert(quaternion_index <= 3); + if (quaternion_index < 0 || quaternion_index > 3) + throw std::invalid_argument("quaternion_index out of range"); } static const unsigned int max_key = 4; @@ -528,7 +532,8 @@ struct Coulombσp_Descr : public Contractable { Coulombσp_Descr() : cartesian_index_(0) {} Coulombσp_Descr(int cartesian_index) : cartesian_index_(cartesian_index) { - assert(cartesian_index >= 0 && cartesian_index <= 2); + if (cartesian_index < 0 || cartesian_index > 2) + throw std::invalid_argument("cartesian_index out of range"); } static const unsigned int max_key = 3; @@ -572,7 +577,8 @@ struct σpσpCoulombσpσp_Descr : public Contractable<σpσpCoulombσpσp_Descr σpσpCoulombσpσp_Descr() : component_index_(0) {} σpσpCoulombσpσp_Descr(int component_index) : component_index_(component_index) { - assert(component_index >= 0 && component_index <= 15); + if (component_index < 0 || component_index > 15) + throw std::invalid_argument("component_index out of range"); } /// 16 components from tensor product of two independent spin spaces: @@ -591,7 +597,7 @@ struct σpσpCoulombσpσp_Descr : public Contractable<σpσpCoulombσpσp_Descr }; // clang-format on const auto qi = component_index(); - if (qi > 15) abort(); + if (qi < 0 || qi > 15) abort(); return std::string("opop_coulomb_opop[") + labels[qi] + "]"; } std::string label() const { return description(); } @@ -642,7 +648,8 @@ struct σpCoulombσp_Descr : public Contractable<σpCoulombσp_Descr> { σpCoulombσp_Descr() : component_index_(0) {} σpCoulombσp_Descr(int component_index) : component_index_(component_index) { - assert(component_index >= 0 && component_index <= 8); + if (component_index < 0 || component_index > 8) + throw std::invalid_argument("component_index out of range"); } static const unsigned int max_key = 9; @@ -654,7 +661,7 @@ struct σpCoulombσp_Descr : public Contractable<σpCoulombσp_Descr> { "symtl_off_xy", "symtl_off_xz", "symtl_off_yz", }; const auto ci = component_index(); - if (ci > 8) abort(); + if (ci < 0 || ci > 8) abort(); return std::string("op_coulomb_op[") + labels[ci] + "]"; } std::string label() const { return description(); } diff --git a/src/bin/libint/rkb_fold_codegen.h b/src/bin/libint/rkb_fold_codegen.h index 43360a91a..28a3b5d7d 100644 --- a/src/bin/libint/rkb_fold_codegen.h +++ b/src/bin/libint/rkb_fold_codegen.h @@ -107,9 +107,9 @@ void rkb_fold_make_children(Factory& factory, const F& a, const F& b, using libint2::prefactor::Scalar; const typename Factory::AuxIndexType zero_m(0u); - std::shared_ptr first_child; - for (std::size_t t = 0; t != terms.size(); ++t) { - const RKBFoldTerm& term = terms[t]; + std::vector> children; + children.reserve(terms.size()); + for (const RKBFoldTerm& term : terms) { F fa{a}; F fb{b}; F fc{c}; @@ -118,28 +118,42 @@ void rkb_fold_make_children(Factory& factory, const F& a, const F& b, if (term.deriv[1] >= 0) fb.deriv().inc(term.deriv[1]); if (term.deriv[2] >= 0) fc.deriv().inc(term.deriv[2]); if (term.deriv[3] >= 0) fd.deriv().inc(term.deriv[3]); - const std::shared_ptr child = - factory.make_child(fa, fb, fc, fd, zero_m); - if (!simple) continue; + children.push_back(factory.make_child(fa, fb, fc, fd, zero_m)); + } + if (!simple) return; - if (t == 0) { - first_child = child; - // a bare child is not an expression: a leading +child is combined with - // the second term below, any other leading term is Scalar(coef) * child - if (term.coef != 1 || terms.size() == 1) expr = Scalar(term.coef) * child; - } else { - const std::shared_ptr lhs = - (t == 1 && terms[0].coef == 1) ? first_child - : std::shared_ptr(expr); - if (term.coef == 1) - expr = lhs + child; - else if (term.coef == -1) - expr = lhs - child; - else - expr = lhs + Scalar(term.coef) * child; + // the order of the children (= the order of the terms) is fixed by + // generate_code(), but the expression may sum them in any order: lead with + // a +1 term, if any, so that no Scalar(-1)*child is needed (a bare child is + // not an expression, hence a lone +1 term is Scalar(1)*child, which is free) + std::size_t lead = 0; + for (std::size_t t = 0; t != terms.size(); ++t) + if (terms[t].coef == 1) { + lead = t; + break; + } + std::shared_ptr lhs = children[lead]; + if (terms[lead].coef != 1) { + expr = Scalar(terms[lead].coef) * children[lead]; + lhs = expr; + nflops += 1; + } else if (terms.size() == 1) { + expr = Scalar(1) * children[lead]; + return; + } + for (std::size_t t = 0; t != terms.size(); ++t) { + if (t == lead) continue; + const auto coef = terms[t].coef; + if (coef == 1) + expr = lhs + children[t]; + else if (coef == -1) + expr = lhs - children[t]; + else { + expr = lhs + Scalar(coef) * children[t]; nflops += 1; } - if (std::abs(term.coef) != 1) nflops += 1; + nflops += 1; + lhs = expr; } } From 91bec9dce8b0fd6e1543c2d2a7cb017ae3ecaec8 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Sun, 25 Oct 2026 23:08:41 -0400 Subject: [PATCH 04/11] Fortran ERI test: compare relative to the magnitude of the shell set The C and Fortran tests compute the primitive data with differently ordered arithmetic; with 2nd-derivative ERIs over d shells (now generated in the CI lane with LIBINT2_ERI_MAX_AM="2;2;2") 2 of ~1400 elements, small differences of large terms, differed by 5e-12 relative to themselves (3e-15 in absolute terms). Measure the difference relative to the largest element of the shell set instead. --- export/fortran/test-eri.cc | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/export/fortran/test-eri.cc b/export/fortran/test-eri.cc index 50cc86306..7cf32dad8 100644 --- a/export/fortran/test-eri.cc +++ b/export/fortran/test-eri.cc @@ -5,6 +5,8 @@ #include #include +#include +#include #include #include "../tests/eri/prep_libint2.h" @@ -498,13 +500,20 @@ bool test_eri_c_f(int &contrdepth, int &am1, double *c1, double *alpha1, const unsigned int n2 = (am2 + 1) * (am2 + 2) / 2; const unsigned int n3 = (am3 + 1) * (am3 + 2) / 2; const unsigned int n4 = (am4 + 1) * (am4 + 2) / 2; + // the C and Fortran codes compute the primitive data with differently + // ordered arithmetic; the resulting roundoff is relative to the magnitude + // of the shell set, not of each element (elements of derivative shell sets + // can be small differences of large terms) + double max_abs = 0.0; + for (unsigned int i = 0; i != n1 * n2 * n3 * n4; ++i) + max_abs = std::max(max_abs, std::abs(eri_shell_set_c[i])); unsigned int nel = 0; for (int a = 0; a < n1; a++) { for (int b = 0; b < n2; b++) { for (int c = 0; c < n3; c++) { for (int d = 0; d < n4; d++, nel++) { const double abs_error = abs(*eri_shell_set_c - *eri_shell_set_f); - const double rel_error = abs_error / abs(*eri_shell_set_c); + const double rel_error = abs_error / max_abs; if (rel_error > threshold) { std::cout << std::setprecision(17) << "Elem " << nel << " di= " << deriv_order From fda4ece485ecb60b9f5bac789c0e43403781cfcf Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Sun, 25 Oct 2026 23:09:19 -0400 Subject: [PATCH 05/11] CHANGES: note PR #450 --- CHANGES | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGES b/CHANGES index e76e390c8..f0127eb7b 100644 --- a/CHANGES +++ b/CHANGES @@ -4,6 +4,10 @@ Following is a brief summary of changes made in each release of Libint. - unreleased: + - PR #450: Engine::compute() throws Engine::lmax_exceeded, in all builds, if a shell exceeds the + max_l of the Engine (previously a Debug-only assert); Engine::lmax_exceeded::lmax_limit() + is the max angular momentum also when thrown by the constructor/set() (previously the + max angular momentum + 1) - PR #399: relativistic restricted kinetic balance (RKB) integrals: 2-body Operator::coulomb_opop, Operator::coulomb_op, Operator::opop_coulomb_opop, Operator::op_coulomb_op (enabled by LIBINT2_ENABLE_RKB_ERI and LIBINT2_ENABLE_RKB_ERI3, off by default; no geometric derivatives yet) From d36ed43baab8eeaef57234cdff6a5b24a91cfde3 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 26 Oct 2026 06:58:23 -0400 Subject: [PATCH 06/11] Engine::compute2: throw lmax_exceeded instead of asserting on the shell AM The per-shell asserts against max_l ran before the new all-builds check, so Debug builds aborted instead of throwing. Replace them with the throw and drop the now-redundant check that followed the canonicalization. --- include/libint2/engine.impl.h | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index 2cd0e4a11..769b51282 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -21,6 +21,7 @@ #ifndef _libint2_src_lib_libint_engineimpl_h_ #define _libint2_src_lib_libint_engineimpl_h_ +#include #include #include "./deriv_map.h" @@ -1347,13 +1348,12 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( tket2.ncontr() == 1) && "generally-contracted shells are not yet supported"); - // angular momentum limit obeyed? can only be fully checked in - // braket-dependent code, here only do a basic test that does not guarantee - // that the shell-set can be computed - assert(tbra1.contr[0].l <= lmax_ && "the angular momentum limit is exceeded"); - assert(tbra2.contr[0].l <= lmax_ && "the angular momentum limit is exceeded"); - assert(tket1.contr[0].l <= lmax_ && "the angular momentum limit is exceeded"); - assert(tket2.contr[0].l <= lmax_ && "the angular momentum limit is exceeded"); + // angular momentum limit obeyed? the scratch, the stack and the core + // evaluator are sized for lmax_, hence check in all builds; the limits of the + // kernels can only be fully checked in braket-dependent code + for (const auto* sh : {&tbra1, &tbra2, &tket1, &tket2}) + if (sh->contr[0].l > lmax_) + throw lmax_exceeded(task_name_, lmax_, sh->contr[0].l); #if LIBINT2_SHELLQUARTET_SET == \ LIBINT2_SHELLQUARTET_SET_STANDARD // standard angular momentum ordering @@ -1437,9 +1437,6 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( // adjust max angular momentum, if needed auto lmax = std::max(std::max(bra1.contr[0].l, bra2.contr[0].l), std::max(ket1.contr[0].l, ket2.contr[0].l)); - // the scratch, the stack and the core evaluator are sized for lmax_, check - // in all builds - if (lmax > lmax_) throw lmax_exceeded(task_name_, lmax_, lmax); const auto lmax_bra = std::max(bra1.contr[0].l, bra2.contr[0].l); const auto lmax_ket = std::max(ket1.contr[0].l, ket2.contr[0].l); From 77d82c4c56a06c4275994015d2e66c0a053153a5 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 26 Oct 2026 07:53:34 -0400 Subject: [PATCH 07/11] Engine: resize primdata_ for the braket rank on set(BraKet)/set(Operator) initialize() resized primdata_ only when given max_nprim, so set(BraKet) and set(Operator) kept the size for the old braket rank; e.g. switching from a 3-center to a 4-center braket overran primdata_ (caught by the new set(BraKet) test, crashing CI). The Engine now records max_nprim_, which initialize() uses, and max_nprim() returns it instead of a value derived from the size of the last computed ShellPair. --- CHANGES | 4 ++++ include/libint2/engine.h | 13 ++++++++----- include/libint2/engine.impl.h | 5 ++++- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/CHANGES b/CHANGES index f0127eb7b..02bfe8005 100644 --- a/CHANGES +++ b/CHANGES @@ -8,6 +8,10 @@ Following is a brief summary of changes made in each release of Libint. max_l of the Engine (previously a Debug-only assert); Engine::lmax_exceeded::lmax_limit() is the max angular momentum also when thrown by the constructor/set() (previously the max angular momentum + 1) + - PR #450: Engine::set(BraKet) and Engine::set(Operator) resize the primitive data for the rank + of the new braket (previously e.g. switching an Engine from a 3-center to a 4-center + braket overran a buffer); Engine::max_nprim() returns the max number of primitives the + Engine was sized for (previously derived from the last computed shell pair) - PR #399: relativistic restricted kinetic balance (RKB) integrals: 2-body Operator::coulomb_opop, Operator::coulomb_op, Operator::opop_coulomb_opop, Operator::op_coulomb_op (enabled by LIBINT2_ENABLE_RKB_ERI and LIBINT2_ENABLE_RKB_ERI3, off by default; no geometric derivatives yet) diff --git a/include/libint2/engine.h b/include/libint2/engine.h index 32e1b12b6..1d3d79cc0 100644 --- a/include/libint2/engine.h +++ b/include/libint2/engine.h @@ -909,6 +909,7 @@ class Engine { primdata_(std::move(other.primdata_)), spbra_(std::move(other.spbra_)), spket_(std::move(other.spket_)), + max_nprim_(other.max_nprim_), stack_size_(other.stack_size_), lmax_(other.lmax_), hard_lmax_(other.hard_lmax_), @@ -941,6 +942,7 @@ class Engine { primdata_(other.primdata_.size()), spbra_(other.spbra_), spket_(other.spket_), + max_nprim_(other.max_nprim_), stack_size_(other.stack_size_), lmax_(other.lmax_), deriv_order_(other.deriv_order_), @@ -964,6 +966,7 @@ class Engine { primdata_ = std::move(other.primdata_); spbra_ = std::move(other.spbra_); spket_ = std::move(other.spket_); + max_nprim_ = other.max_nprim_; stack_size_ = other.stack_size_; lmax_ = other.lmax_; hard_lmax_ = other.hard_lmax_; @@ -997,6 +1000,7 @@ class Engine { primdata_.resize(other.primdata_.size()); spbra_ = other.spbra_; spket_ = other.spket_; + max_nprim_ = other.max_nprim_; stack_size_ = other.stack_size_; lmax_ = other.lmax_; deriv_order_ = other.deriv_order_; @@ -1075,17 +1079,14 @@ class Engine { } /// @return the maximum number of primitives that this engine can handle - std::size_t max_nprim() const { - assert(spbra_.primpairs.size() == spket_.primpairs.size()); - return static_cast(std::sqrt(spbra_.primpairs.size())); - } + std::size_t max_nprim() const { return max_nprim_; } /// reset the maximum number of primitives /// @param[in] n the maximum number of primitives /// @note left unchanged if the value returned by Engine::max_nprim is greater /// than @c n Engine& set_max_nprim(std::size_t n) { - if (n * n > spbra_.primpairs.size()) { + if (n > max_nprim_) { spbra_.resize(n); spket_.resize(n); initialize(n); @@ -1338,6 +1339,8 @@ class Engine { BraKet braket_; std::vector primdata_; ShellPair spbra_, spket_; + // primdata_ is sized for max_nprim_ primitives per shell; set by initialize() + std::size_t max_nprim_ = 0; size_t stack_size_; // amount allocated by libint2_init_xxx in // primdata_[0].stack int lmax_; diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index 769b51282..17dcc5dfb 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -721,7 +721,10 @@ __libint2_engine_inline void Engine::initialize(size_t max_nprim) { // initialize braket, if needed if (braket_ == BraKet::invalid) braket_ = default_braket(oper_); - if (max_nprim != 0) primdata_.resize(std::pow(max_nprim, braket_rank())); + // the size of primdata_ depends on the rank of the braket, hence resize it + // also when called without max_nprim, e.g. by set(BraKet) or set(Operator) + if (max_nprim != 0) max_nprim_ = max_nprim; + primdata_.resize(std::pow(max_nprim_, braket_rank())); // initialize targets { From eaf4fe993d2bb8d74f12e9c4667b63ffafd3fa0f Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 26 Oct 2026 09:29:28 -0400 Subject: [PATCH 08/11] Engine::compute2: evaluate the derivative order at compile time intrinsic_deriv_order() dispatches on the runtime operator, and compute2 called it in the loop over primitives; the operator and the derivative order are template parameters of compute2, hence use a constexpr. --- include/libint2/engine.impl.h | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index 17dcc5dfb..1b254eaaa 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -1340,6 +1340,10 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( (tspbra != nullptr && tspket != nullptr)) && "Engine::compute2 -- expects zero or two ShellPair objects"); assert(screening_method_ != ScreeningMethod::Invalid); + // the derivative order of the integrals over the operator, known at compile + // time, hence evaluated once rather than per primitive + constexpr auto total_deriv_order = + der + operator_traits::intrinsic_deriv_order; if (buildfnptrs_ == nullptr) throw missing_kernel(); // @@ -1599,7 +1603,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( const scalar_type rho = gammap * gammaq * oogammapq; const scalar_type T = PQ2 * rho; auto* gm_ptr = &(primdata.LIBINT_T_SS_EREP_SS(0)[0]); - const auto mmax = l + deriv_order_ + intrinsic_deriv_order(); + const auto mmax = l + total_deriv_order; if (!skip_core_ints) { // Simple core-eval dispatch: just `core_eval_ptr->eval(gm_ptr, T, mmax)`. @@ -1849,7 +1853,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( const auto Wz = (gammap_o_gammapgammaq * P[2] + gammaq_o_gammapgammaq * Q[2]); - if (deriv_order_ + intrinsic_deriv_order() > 0 || lmax_bra > 0) { + if (total_deriv_order > 0 || lmax_bra > 0) { #if LIBINT2_DEFINED(eri, WP_x) primdata.WP_x[0] = Wx - P[0]; #endif @@ -1860,7 +1864,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( primdata.WP_z[0] = Wz - P[2]; #endif } - if (deriv_order_ + intrinsic_deriv_order() > 0 || lmax_ket > 0) { + if (total_deriv_order > 0 || lmax_ket > 0) { #if LIBINT2_DEFINED(eri, WQ_x) primdata.WQ_x[0] = Wx - Q[0]; #endif @@ -1944,7 +1948,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( #endif // prefactors for derivative ERI relations - if (deriv_order_ + intrinsic_deriv_order() > 0) { + if (total_deriv_order > 0) { #if LIBINT2_DEFINED(eri, alpha1_rho_over_zeta2) primdata.alpha1_rho_over_zeta2[0] = alpha0 * (oogammap * gammaq_o_gammapgammaq); @@ -2027,8 +2031,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( } // compute directly (ss|ss) - const auto compute_directly = - lmax == 0 && deriv_order_ == 0 && intrinsic_deriv_order() == 0; + const auto compute_directly = lmax == 0 && total_deriv_order == 0; if (compute_directly) { #ifdef LIBINT2_ENGINE_TIMERS From dcd020d975bbb7befa9d7e6c1a2441e2df7645ce Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 26 Oct 2026 09:29:29 -0400 Subject: [PATCH 09/11] Compiler: a recurrence relation whose target is its child is an alias MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit RecurrenceRelation::is_alias() (false by default) declares that the target is identical, element by element and in the same layout, to the only child. handle_trivial_nodes() then makes the vertex refer to the child, provided the child has no other parent and is not a target, so no code is generated for the RR; update_func_names() and ExtractRR skip referring vertices, so neither is the RR's function included or generated. CR_11_Coulombσp_11 (one derivative integral with unit coefficient per component) is an alias: the 19 3coulomb_op prereq kernels of the RKB test configuration accumulate the derivative integrals directly instead of copying them first, and CR_Coulombop_{0,1,2} are no longer generated. --- .../bin/libint/comp_11_Coulomb\317\203p_11.h" | 11 +++++-- src/bin/libint/dg.cc | 29 +++++++++++++++++++ src/bin/libint/extract.cc | 4 +++ src/bin/libint/rr.h | 7 +++++ 4 files changed, 49 insertions(+), 2 deletions(-) diff --git "a/src/bin/libint/comp_11_Coulomb\317\203p_11.h" "b/src/bin/libint/comp_11_Coulomb\317\203p_11.h" index fd85b01f4..a0450eb36 100644 --- "a/src/bin/libint/comp_11_Coulomb\317\203p_11.h" +++ "b/src/bin/libint/comp_11_Coulomb\317\203p_11.h" @@ -25,11 +25,10 @@ #include #include #include +#include #include #include -#include - #include namespace libint2 { @@ -64,6 +63,14 @@ class CR_11_Coulombσp_11 : public GenericRecurrenceRelation< static bool directional() { return false; } + /// each Cartesian component is a single derivative integral with unit + /// coefficient, hence the target is the child (see fold_terms()) + bool is_alias() const override { + if (this->num_children() != 1) return false; + const auto terms = fold_terms(target_->oper()->descr().cartesian_index()); + return terms.size() == 1 && terms[0].coef == 1; + } + private: using ParentType::is_simple; using ParentType::target_; diff --git a/src/bin/libint/dg.cc b/src/bin/libint/dg.cc index 0955cbb05..3d9b901ac 100644 --- a/src/bin/libint/dg.cc +++ b/src/bin/libint/dg.cc @@ -846,6 +846,32 @@ void DirectedGraph::handle_trivial_nodes( } } + // Is the exit arc an RR whose target is identical to its child (e.g. + // CR_11_Coulombσp_11)? Then make the vertex refer to the child, so that the + // child is computed in its place and no code is generated for the RR. + // N.B. the child must have no other parent: a referring vertex is never + // allocated, hence never tags its child, and freeing it (once its parents + // are allocated) frees the child; if it is a target, the child takes over + // its address. + if (!(vptr)->refers_to_another()) { + std::shared_ptr arc_cast = + std::dynamic_pointer_cast(arc); + if (arc_cast && arc_cast->rr()->is_alias()) { + std::shared_ptr child = arc->dest(); + // the vertex and the child must not be pinned to different locations + const bool location_conflict = + ((vptr)->address_set() || (vptr)->symbol_set()) && + (child->address_set() || child->symbol_set()); + if (!location_conflict && !child->is_a_target() && + child->num_entry_arcs() == 1) { + if ((vptr)->address_set()) child->set_address((vptr)->address()); + vptr->refer_this_to(child); + vptr->reset_symbol(); + continue; + } + } + } + // Is the exit arc DGArcDirect? { std::shared_ptr arc_cast = @@ -2270,6 +2296,9 @@ void DirectedGraph::update_func_names() { typedef vertices::iterator iter; for (iter v = stack_.begin(); v != stack_.end(); ++v) { const ver_ptr& vptr = vertex_ptr(*v); + // no code is generated for a vertex that refers to another (see + // handle_trivial_nodes()), hence no function is called for it + if ((vptr)->refers_to_another()) continue; // for every vertex with children if ((vptr)->num_exit_arcs() > 0) { // if it must be computed using a RR diff --git a/src/bin/libint/extract.cc b/src/bin/libint/extract.cc index 24b86252f..c395a5a1f 100644 --- a/src/bin/libint/extract.cc +++ b/src/bin/libint/extract.cc @@ -79,6 +79,10 @@ const ExtractExternSymbols::Symbols& ExtractExternSymbols::symbols() const { //// void ExtractRR::operator()(const VertexPtr& v) { + // no code is generated for a vertex that refers to another (e.g. a vertex + // whose RR is an alias, see DirectedGraph::handle_trivial_nodes()), hence + // neither is its RR used + if (v->refers_to_another()) return; if (v->num_exit_arcs() != 0) { std::shared_ptr arc = *(v->first_exit_arc()); std::shared_ptr arc_rr = diff --git a/src/bin/libint/rr.h b/src/bin/libint/rr.h index 1f46b65c1..bacc8172f 100644 --- a/src/bin/libint/rr.h +++ b/src/bin/libint/rr.h @@ -125,6 +125,13 @@ class RecurrenceRelation inlined and optimized. */ virtual bool is_simple() const = 0; + /** + Returns true if the target is identical, element by element and in the + same layout, to the only child; then the vertex of the target is made to + refer to the child (see DirectedGraph::handle_trivial_nodes()), hence no + code is generated for this recurrence relation. + */ + virtual bool is_alias() const { return false; } /** Returns true is the type of target and all children are exactly the same */ From 5fb646929dbe2e7bcaf3657afded9494f62316bc Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 26 Oct 2026 09:29:30 -0400 Subject: [PATCH 10/11] Compiler test: DerivGaussV2; fail on an exception MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit test9 checks the children of direct differentiation and of translational invariance, that the generated code is shared by the spectators but not by different shells or directions, builds shell sets that call both kinds of DerivGaussV2 functions, and that CR_11_Coulombσp_11 is an alias. The test driver returned 0 even if a test threw, hence could never fail. --- src/bin/libint/test.cc | 149 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 149 insertions(+) diff --git a/src/bin/libint/test.cc b/src/bin/libint/test.cc index 9f9d0c061..ca47e3d0d 100644 --- a/src/bin/libint/test.cc +++ b/src/bin/libint/test.cc @@ -53,6 +53,7 @@ void test5(); void test6(); void test7(); void test8(); +void test9(); void test_cgshell_iter(const CGShell& sh); template @@ -88,6 +89,7 @@ int main(int argc, char* argv[]) { << " WARNING! Caught a standard exception:" << endl << " " << a.what() << endl << endl; + return 1; } return 0; } @@ -148,6 +150,9 @@ int try_main(int argc, char* argv[]) { #if 1 RunTest(test8, "contracted G12 integral build"); #endif +#if 1 + RunTest(test9, "DerivGaussV2 relation"); +#endif return 0; } @@ -517,4 +522,148 @@ void test8() { cparams->contracted_targets(contracted_targets_old_value); } +// DerivGaussV2: the children of direct differentiation and of translational +// invariance, and the labels by which the RRStack shares the generated code +void test9() { + typedef TwoPRep_11_11_sq IType; + auto check = [](bool ok, const std::string& what) { + if (!ok) throw ProgrammingError("DerivGaussV2 test: " + what); + }; + // @return whether child @p i of @p rr is @p iset + auto child_is = [](const std::shared_ptr& rr, + unsigned int i, const std::shared_ptr& iset) { + return rr->rr_child(i) == std::static_pointer_cast(iset); + }; + // @return shell of angular momentum @p l (primitive, unless @p contracted) + // with derivative quanta @p dx, @p dy, @p dz + auto shell = [](unsigned int l, unsigned int dx = 0, unsigned int dy = 0, + unsigned int dz = 0, bool contracted = false) { + CGShell sh(l); + if (contracted) + sh.contract(); + else + sh.uncontract(); + if (dx) sh.deriv().inc(0, dx); + if (dy) sh.deriv().inc(1, dy); + if (dz) sh.deriv().inc(2, dz); + return sh; + }; + const mType m0(0u); + + // direct differentiation of a (particle 0, bra), in x: + // ∂_Ax (p_x s|s s) = 2α (d s|s s) - 1 (s s|s s) + { + auto target = + IType::Instance(shell(1, 1), shell(0), shell(0), shell(0), m0); + auto rr = DerivV2_a_11_TwoPRep_11_sh::Instance(target, 0); + check(rr != nullptr, "no relation for (p_x s|s s)"); + check(rr->num_children() == 2, "(p_x s|s s) must have 2 children"); + check(child_is(rr, 0, + IType::Instance(shell(2), shell(0), shell(0), shell(0), m0)), + "child 0 of (p_x s|s s) must be (d s|s s)"); + check(child_is(rr, 1, + IType::Instance(shell(0), shell(0), shell(0), shell(0), m0)), + "child 1 of (p_x s|s s) must be (s s|s s)"); + // s has no a-1 term + auto target_s = + IType::Instance(shell(0, 1), shell(0), shell(0), shell(0), m0); + auto rr_s = DerivV2_a_11_TwoPRep_11_sh::Instance(target_s, 0); + check(rr_s != nullptr && rr_s->num_children() == 1, + "(s_x s|s s) must have 1 child"); + check(child_is(rr_s, 0, + IType::Instance(shell(1), shell(0), shell(0), shell(0), m0)), + "child 0 of (s_x s|s s) must be (p s|s s)"); + // no derivative quanta in y, no relation + check(DerivV2_a_11_TwoPRep_11_sh::Instance(target, 1) == nullptr, + "(p_x s|s s) must not be differentiated in y"); + // direct differentiation applies to primitives only + auto target_c = + IType::Instance(shell(1, 1, 0, 0, true), shell(0, 0, 0, 0, true), + shell(0, 0, 0, 0, true), shell(0, 0, 0, 0, true), m0); + check(DerivV2_a_11_TwoPRep_11_sh::Instance(target_c, 0) == nullptr, + "contracted (p_x s|s s) must not be differentiated directly"); + } + + // code sharing: the label depends on the differentiated shell and the + // direction only, not on the spectators + { + auto rr1 = DerivV2_a_11_TwoPRep_11_sh::Instance( + IType::Instance(shell(1, 1), shell(0), shell(0), shell(0), m0), 0); + auto rr2 = DerivV2_a_11_TwoPRep_11_sh::Instance( + IType::Instance(shell(1, 1), shell(2), shell(1), shell(1), m0), 0); + auto rr3 = DerivV2_a_11_TwoPRep_11_sh::Instance( + IType::Instance(shell(2, 1), shell(0), shell(0), shell(0), m0), 0); + check(rr1 != nullptr && rr2 != nullptr && rr3 != nullptr, + "no relation for (p_x s|s s), (p_x d|p p) or (d_x s|s s)"); + check(rr1->label() == rr2->label(), + "(p_x s|s s) and (p_x d|p p) must share the code: " + rr1->label() + + " vs " + rr2->label()); + check(rr1->label() != rr3->label(), + "(p_x s|s s) and (d_x s|s s) must not share the code"); + // a shell with derivative quanta in x and y is differentiated in either + // direction, with different children: the labels must differ + auto target_xy = + IType::Instance(shell(1, 1, 1), shell(0), shell(0), shell(0), m0); + auto rr_x = DerivV2_a_11_TwoPRep_11_sh::Instance(target_xy, 0); + auto rr_y = DerivV2_a_11_TwoPRep_11_sh::Instance(target_xy, 1); + check(rr_x != nullptr && rr_y != nullptr, + "(p_xy s|s s) must be differentiated in x and y"); + check(rr_x->label() != rr_y->label(), + "the x and y derivatives of (p_xy s|s s) must not share the code"); + } + + // translational invariance (particle 1, bra = c), also for contracted + // shells: ∂_Cx (a b|c d) = -∂_Ax (a b|c d) - ∂_Bx (a b|c d) - ∂_Dx (a b|c d) + { + typedef DerivGaussV2 TIType; + const bool c = true; + auto sh = [&](unsigned int l, unsigned int dx = 0) { + return shell(l, dx, 0, 0, c); + }; + auto target = IType::Instance(sh(1), sh(0), sh(1, 1), sh(1), m0); + auto rr = TIType::Instance(target, 0); + check(rr != nullptr, "no translational invariance relation"); + check(rr->num_children() == 3, + "translational invariance must have 3 children"); + check(child_is(rr, 0, IType::Instance(sh(1, 1), sh(0), sh(1), sh(1), m0)), + "child 0 must be (p_x s|p p)"); + check(child_is(rr, 1, IType::Instance(sh(1), sh(0, 1), sh(1), sh(1), m0)), + "child 1 must be (p s_x|p p)"); + check(child_is(rr, 2, IType::Instance(sh(1), sh(0), sh(1), sh(1, 1), m0)), + "child 2 must be (p s|p p_x)"); + check(rr->label().rfind("DerivGaussV2TI", 0) == 0, + "translational invariance label must start with DerivGaussV2TI: " + + rr->label()); + } + + // build contracted shell sets, without unrolling, hence with calls to the + // generated DerivGaussV2 functions: differentiation of a (direct, shared by + // the spectators) and of c (translational invariance) + { + CGShell::set_contracted_default_value(true); + const bool contracted_targets_old_value = cparams->contracted_targets(); + cparams->contracted_targets(true); + const bool c = true; + RunBuildTest(shell(1, 1, 0, 0, c), shell(1, 0, 0, 0, c), + shell(0, 0, 0, 0, c), shell(1, 0, 0, 0, c), 0, + use_quartets); + RunBuildTest(shell(1, 0, 0, 0, c), shell(0, 0, 0, 0, c), + shell(1, 0, 1, 0, c), shell(0, 0, 0, 0, c), 0, + use_quartets); + cparams->contracted_targets(contracted_targets_old_value); + } + + // CR_11_Coulombσp_11 is an alias (target = child, see + // DirectedGraph::handle_trivial_nodes()), CR_11_Coulombσpσp_11 is not + { + for (int k = 0; k != 3; ++k) { + auto target = Coulombσp_11_11_sq::Instance( + shell(1), shell(0), shell(1), shell(1), m0, Coulombσp_Descr(k)); + auto rr = CR_11_Coulombσp_11_sh::Instance(target, 0); + check(rr != nullptr && rr->num_children() == 1 && rr->is_alias(), + "CR_11_Coulombσp_11 must be an alias"); + } + } +} + }; // namespace From 54c8e369f6b49c6873afb0ef2a6caeb9fb4f4a46 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Mon, 26 Oct 2026 09:29:31 -0400 Subject: [PATCH 11/11] CHANGES: coulomb_op kernels compute in place --- CHANGES | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGES b/CHANGES index 02bfe8005..c48596a6f 100644 --- a/CHANGES +++ b/CHANGES @@ -12,6 +12,9 @@ Following is a brief summary of changes made in each release of Libint. of the new braket (previously e.g. switching an Engine from a 3-center to a 4-center braket overran a buffer); Engine::max_nprim() returns the max number of primitives the Engine was sized for (previously derived from the last computed shell pair) + - PR #450: the generated Operator::coulomb_op kernels compute each component in place, without + an element-wise copy of the derivative integrals (the compiler makes the target of a + recurrence relation that is an alias of its child refer to the child) - PR #399: relativistic restricted kinetic balance (RKB) integrals: 2-body Operator::coulomb_opop, Operator::coulomb_op, Operator::opop_coulomb_opop, Operator::op_coulomb_op (enabled by LIBINT2_ENABLE_RKB_ERI and LIBINT2_ENABLE_RKB_ERI3, off by default; no geometric derivatives yet)