From 2e2aae4a8380a6a99b4b35499f110d152a754535 Mon Sep 17 00:00:00 2001 From: Andreas Koehn Date: Fri, 2 Oct 2026 17:17:41 +0200 Subject: [PATCH 1/4] Add Operator::delcgtga_delcgtgb to Engine Two-geminal generalization of delcgtg2: (Delta . CGTG a) . (Delta . CGTG b), with params std::array{a, b}. Co-Authored-By: Claude Opus 5.5 --- include/libint2/engine.h | 13 ++++++++++++- include/libint2/engine.impl.h | 35 ++++++++++++++++++++++++++++++++++- 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/include/libint2/engine.h b/include/libint2/engine.h index f07e8096c..1ce0f189d 100644 --- a/include/libint2/engine.h +++ b/include/libint2/engine.h @@ -225,6 +225,10 @@ enum class Operator { stg_x_coulomb, /// alias for Operator::stg_x_coulomb yukawa = stg_x_coulomb, + /// (Delta . contracted Gaussian geminal a) . (Delta . contracted Gaussian + /// geminal b), i.e. the two-geminal generalization of delcgtg2; params are + /// std::array{a, b} + delcgtga_delcgtgb, // do not modify this invalid = -1, // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!keep this @@ -232,7 +236,7 @@ enum class Operator { first_1body_oper = overlap, last_1body_oper = op_q_gau_op, first_2body_oper = delta, - last_2body_oper = stg_x_coulomb, + last_2body_oper = delcgtga_delcgtgb, first_oper = first_1body_oper, last_oper = last_2body_oper }; @@ -451,6 +455,13 @@ struct operator_traits template <> struct operator_traits : public detail::cgtg_operator_traits<2> {}; +template <> +struct operator_traits + : public detail::cgtg_operator_traits<2> { + /// the two geminals a and b + typedef std::array oper_params_type; + static oper_params_type default_params() { return oper_params_type{}; } +}; template <> struct operator_traits diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index f2f7cedcd..c94c92838 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -95,7 +95,8 @@ typename std::remove_all_extents::type* to_ptr1(T (&a)[N]) { (eri, /* erfx_coulomb */ \ (eri, /* stg */ \ (eri, /* yukawa */ \ - BOOST_PP_NIL)))))))))))))))))))))))) + (eri, /* delcgtga_delcgtgb */ \ + BOOST_PP_NIL))))))))))))))))))))))))) #define BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE \ BOOST_PP_MAKE_TUPLE(BOOST_PP_LIST_SIZE(BOOST_PP_NBODY_OPERATOR_LIST)) @@ -942,6 +943,25 @@ __libint2_engine_inline void Engine::init_core_ints_params(const any& params) { std::to_string(std::get<1>(q_gau_params).size()) + ")"); } core_ints_params_ = params; + } else if (oper_ == Operator::delcgtga_delcgtgb) { + // (\Del g_a) \cdot (\Del g_b) with g_a = \sum_i c_i exp(-a_i r_12^2), + // g_b = \sum_j d_j exp(-b_j r_12^2): + // = \sum_ij 4 a_i b_j c_i d_j r_{12}^2 exp(-(a_i+b_j) r_{12}^2) + // i.e. product of the two geminals, each coefficient scaled by 4 a_i b_j + const auto& oparams = any_cast::oper_params_type&>(params); + const auto& ga = oparams[0]; + const auto& gb = oparams[1]; + operator_traits::oper_params_type core_ints_params; + core_ints_params.reserve(ga.size() * gb.size()); + for (const auto& pa : ga) + for (const auto& pb : gb) { + const auto gexp = pa.first + pb.first; + const auto gcoeff_rescaled = + 4 * pa.first * pb.first * pa.second * pb.second; + core_ints_params.push_back(std::make_pair(gexp, gcoeff_rescaled)); + } + core_ints_params_ = core_ints_params; } else { core_ints_params_ = params; } @@ -1509,6 +1529,19 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( core_ints_params_); core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params); } break; + case Operator::delcgtga_delcgtgb: { + // core_ints_params_ holds the product geminal, see + // init_core_ints_params + const auto& core_eval_ptr = + any_cast&>(core_eval_pack_) + .first(); + const auto& core_ints_params = + any_cast::oper_params_type&>( + core_ints_params_); + core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params); + } break; case Operator::delta: { const auto& core_eval_ptr = any_cast< From 2354c2a6fe80f2e283eb0d346a1ca7a9b5b7fac7 Mon Sep 17 00:00:00 2001 From: Andreas Koehn Date: Sat, 24 Oct 2026 23:45:44 +0200 Subject: [PATCH 2/4] Add unit test for Operator::delcgtga_delcgtgb Checks against delcgtg2: for a = b the integrals must reproduce delcgtg2(a); for a != b they must be symmetric in (a, b) and satisfy the polarization identity (Del g_a).(Del g_b) = (delcgtg2(g_a + g_b) - delcgtg2(g_a) - delcgtg2(g_b)) / 2 Co-Authored-By: Claude Opus 5.5 --- export/tests/unit/test-2body.cc | 85 +++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) diff --git a/export/tests/unit/test-2body.cc b/export/tests/unit/test-2body.cc index fd602a910..9a7a471ed 100644 --- a/export/tests/unit/test-2body.cc +++ b/export/tests/unit/test-2body.cc @@ -122,6 +122,91 @@ TEST_CASE("Slater/Yukawa integrals", "[engine][2-body]") { } } +TEST_CASE("delcgtga_delcgtgb integrals", "[engine][2-body]") { + std::vector obs{ + // pseudorandom s + Shell{{1.5, 0.4}, {{0, true, {1.0, 0.5}}}, {{0.0, 0.0, 0.0}}}, + // pseudorandom p + Shell{{1.0, 3.0}, {{1, true, {1.0, 0.3}}}, {{0.5, -0.3, 0.2}}}, + // pseudorandom d + Shell{{2.0, 5.0}, {{2, true, {1.0, 0.2}}}, {{1.0, 1.0, 1.0}}}, + }; + const auto max_nprim = libint2::max_nprim(obs); + const auto max_l = libint2::max_l(obs); + + // two unrelated pseudorandom geminals + using geminal_t = std::vector>; + const geminal_t ga = {{0.3, 0.4}, {1.2, 0.25}, {4.5, -0.1}}; + const geminal_t gb = {{0.7, 0.6}, {2.5, -0.2}}; + geminal_t ga_plus_gb = ga; + ga_plus_gb.insert(ga_plus_gb.end(), gb.begin(), gb.end()); + + if (LIBINT2_MAX_AM_eri >= max_l) { + auto engine = Engine(Operator::delcgtga_delcgtgb, max_nprim, max_l); + // N.B. delcgtg2 has no default params, must construct with params + auto engine_ref = Engine(Operator::delcgtg2, max_nprim, max_l, 0, + std::numeric_limits::epsilon(), ga); + + // computes all shell quartets of obs with engine e + auto compute_all = [&](Engine& e) { + std::vector ints; + const auto nshell = obs.size(); + for (int s1 = 0; s1 != nshell; ++s1) { + for (int s2 = 0; s2 != nshell; ++s2) { + for (int s3 = 0; s3 != nshell; ++s3) { + for (int s4 = 0; s4 != nshell; ++s4) { + e.compute(obs[s1], obs[s2], obs[s3], obs[s4]); + const auto setsize = obs[s1].size() * obs[s2].size() * + obs[s3].size() * obs[s4].size(); + if (e.results()[0] == nullptr) + ints.insert(ints.end(), setsize, 0.0); + else + ints.insert(ints.end(), e.results()[0], + e.results()[0] + setsize); + } + } + } + } + return ints; + }; + + SECTION("a = b reproduces delcgtg2") { + engine.set_params(std::array{ga, ga}); + engine_ref.set_params(ga); + const auto ints = compute_all(engine); + const auto ints_ref = compute_all(engine_ref); + REQUIRE(ints.size() == ints_ref.size()); + for (std::size_t i = 0; i != ints.size(); ++i) + REQUIRE(ints[i] == Approx(ints_ref[i]).margin(1e-12)); + } + + SECTION("a != b: symmetric and satisfies polarization identity") { + // (Del g_a).(Del g_b) = ((Del (g_a + g_b))^2 - (Del g_a)^2 - + // (Del g_b)^2) / 2 + engine.set_params(std::array{ga, gb}); + const auto ints_ab = compute_all(engine); + engine.set_params(std::array{gb, ga}); + const auto ints_ba = compute_all(engine); + engine_ref.set_params(ga_plus_gb); + const auto ints_apb = compute_all(engine_ref); + engine_ref.set_params(ga); + const auto ints_a = compute_all(engine_ref); + engine_ref.set_params(gb); + const auto ints_b = compute_all(engine_ref); + + bool nontrivial = false; + for (std::size_t i = 0; i != ints_ab.size(); ++i) { + REQUIRE(ints_ab[i] == Approx(ints_ba[i]).margin(1e-12)); + REQUIRE(ints_ab[i] == + Approx((ints_apb[i] - ints_a[i] - ints_b[i]) / 2) + .margin(1e-10)); + if (std::abs(ints_ab[i]) > 1e-3) nontrivial = true; + } + REQUIRE(nontrivial); + } + } +} + // see https://github.com/evaleev/libint/issues/216 TEST_CASE_METHOD(libint2::unit::DefaultFixture, "bra-ket permutation", "[engine][2-body]") { From cff28f8427a4372eda056c45ffc2418f371194a7 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Sun, 25 Oct 2026 10:47:52 -0400 Subject: [PATCH 3/4] delcgtga_delcgtgb: require params, validate them, fix stale params after set(Operator) - drop the empty default params of delcgtga_delcgtgb (like delcgtg2), so an engine constructed without params throws instead of silently producing zeros - init_core_ints_params: descriptive errors for missing/wrong-typed params and negative geminal exponents for delcgtg2 and delcgtga_delcgtgb - Engine::set(Operator) resets core_ints_params_, so params of the previous operator (e.g. delcgtg2's geminal) are no longer used by the new operator - share the compute2 case of delcgtg2 and delcgtga_delcgtgb - fix the doc of delcgtga_delcgtgb (gradient, not Laplacian; relation to delcgtg2) - test: 1st derivatives, required/validated params, set(Operator) --- export/tests/unit/test-2body.cc | 118 +++++++++++++++++++++++--------- include/libint2/engine.h | 15 ++-- include/libint2/engine.impl.h | 113 +++++++++++++++++------------- 3 files changed, 162 insertions(+), 84 deletions(-) diff --git a/export/tests/unit/test-2body.cc b/export/tests/unit/test-2body.cc index 9a7a471ed..d25bcc243 100644 --- a/export/tests/unit/test-2body.cc +++ b/export/tests/unit/test-2body.cc @@ -141,53 +141,73 @@ TEST_CASE("delcgtga_delcgtgb integrals", "[engine][2-body]") { geminal_t ga_plus_gb = ga; ga_plus_gb.insert(ga_plus_gb.end(), gb.begin(), gb.end()); - if (LIBINT2_MAX_AM_eri >= max_l) { - auto engine = Engine(Operator::delcgtga_delcgtgb, max_nprim, max_l); - // N.B. delcgtg2 has no default params, must construct with params - auto engine_ref = Engine(Operator::delcgtg2, max_nprim, max_l, 0, - std::numeric_limits::epsilon(), ga); - - // computes all shell quartets of obs with engine e - auto compute_all = [&](Engine& e) { - std::vector ints; - const auto nshell = obs.size(); - for (int s1 = 0; s1 != nshell; ++s1) { - for (int s2 = 0; s2 != nshell; ++s2) { - for (int s3 = 0; s3 != nshell; ++s3) { - for (int s4 = 0; s4 != nshell; ++s4) { - e.compute(obs[s1], obs[s2], obs[s3], obs[s4]); - const auto setsize = obs[s1].size() * obs[s2].size() * - obs[s3].size() * obs[s4].size(); - if (e.results()[0] == nullptr) + if (LIBINT2_MAX_AM_eri < max_l) { + WARN("delcgtga_delcgtgb integrals: skipped, LIBINT2_MAX_AM_eri < " + << max_l); + return; + } + // also test 1st derivatives, if available + const int max_deriv_order = +#if LIBINT2_DERIV_ERI_ORDER > 0 + (LIBINT2_MAX_AM_eri1 >= max_l) ? 1 : 0; +#else + 0; +#endif + const auto eps = std::numeric_limits::epsilon(); + using params_t = std::array; + + // computes all shell quartets of obs (all targets, i.e. all derivatives) + // with engine e + auto compute_all = [&](Engine &e) { + std::vector ints; + const auto nshell = obs.size(); + for (std::size_t s1 = 0; s1 != nshell; ++s1) { + for (std::size_t s2 = 0; s2 != nshell; ++s2) { + for (std::size_t s3 = 0; s3 != nshell; ++s3) { + for (std::size_t s4 = 0; s4 != nshell; ++s4) { + e.compute(obs[s1], obs[s2], obs[s3], obs[s4]); + const auto setsize = obs[s1].size() * obs[s2].size() * + obs[s3].size() * obs[s4].size(); + for (const auto *target : e.results()) { + if (target == nullptr) ints.insert(ints.end(), setsize, 0.0); else - ints.insert(ints.end(), e.results()[0], - e.results()[0] + setsize); + ints.insert(ints.end(), target, target + setsize); } } } } - return ints; - }; + } + return ints; + }; - SECTION("a = b reproduces delcgtg2") { - engine.set_params(std::array{ga, ga}); - engine_ref.set_params(ga); + SECTION("a = b reproduces delcgtg2") { + for (int d = 0; d <= max_deriv_order; ++d) { + CAPTURE(d); + auto engine = Engine(Operator::delcgtga_delcgtgb, max_nprim, max_l, d, + eps, params_t{ga, ga}); + auto engine_ref = + Engine(Operator::delcgtg2, max_nprim, max_l, d, eps, ga); const auto ints = compute_all(engine); const auto ints_ref = compute_all(engine_ref); REQUIRE(ints.size() == ints_ref.size()); for (std::size_t i = 0; i != ints.size(); ++i) REQUIRE(ints[i] == Approx(ints_ref[i]).margin(1e-12)); } + } - SECTION("a != b: symmetric and satisfies polarization identity") { - // (Del g_a).(Del g_b) = ((Del (g_a + g_b))^2 - (Del g_a)^2 - - // (Del g_b)^2) / 2 - engine.set_params(std::array{ga, gb}); + SECTION("a != b: symmetric and satisfies polarization identity") { + // (Del g_a).(Del g_b) = ((Del (g_a + g_b))^2 - (Del g_a)^2 - + // (Del g_b)^2) / 2 + for (int d = 0; d <= max_deriv_order; ++d) { + CAPTURE(d); + auto engine = Engine(Operator::delcgtga_delcgtgb, max_nprim, max_l, d, + eps, params_t{ga, gb}); + auto engine_ref = + Engine(Operator::delcgtg2, max_nprim, max_l, d, eps, ga_plus_gb); const auto ints_ab = compute_all(engine); - engine.set_params(std::array{gb, ga}); + engine.set_params(params_t{gb, ga}); const auto ints_ba = compute_all(engine); - engine_ref.set_params(ga_plus_gb); const auto ints_apb = compute_all(engine_ref); engine_ref.set_params(ga); const auto ints_a = compute_all(engine_ref); @@ -197,14 +217,44 @@ TEST_CASE("delcgtga_delcgtgb integrals", "[engine][2-body]") { bool nontrivial = false; for (std::size_t i = 0; i != ints_ab.size(); ++i) { REQUIRE(ints_ab[i] == Approx(ints_ba[i]).margin(1e-12)); - REQUIRE(ints_ab[i] == - Approx((ints_apb[i] - ints_a[i] - ints_b[i]) / 2) - .margin(1e-10)); + REQUIRE( + ints_ab[i] == + Approx((ints_apb[i] - ints_a[i] - ints_b[i]) / 2).margin(1e-10)); if (std::abs(ints_ab[i]) > 1e-3) nontrivial = true; } REQUIRE(nontrivial); } } + + SECTION("params are required and validated") { + // no default params + REQUIRE_THROWS_AS(Engine(Operator::delcgtga_delcgtgb, max_nprim, max_l), + std::logic_error); + // a single geminal is not valid + REQUIRE_THROWS( + Engine(Operator::delcgtga_delcgtgb, max_nprim, max_l, 0, eps, ga)); + auto engine = Engine(Operator::delcgtga_delcgtgb, max_nprim, max_l, 0, eps, + params_t{ga, gb}); + REQUIRE_THROWS_AS(engine.set_params(ga), std::logic_error); + // negative exponents are rejected + const geminal_t gneg = {{-0.5, 1.0}}; + REQUIRE_THROWS_AS(engine.set_params(params_t{ga, gneg}), std::logic_error); + } + + SECTION("set(Operator) does not reuse params of the previous operator") { + auto engine = Engine(Operator::delcgtg2, max_nprim, max_l, 0, eps, ga); + engine.set(Operator::delcgtga_delcgtgb); + // delcgtg2's params must not be used, hence no params available + REQUIRE_THROWS(compute_all(engine)); + engine.set_params(params_t{ga, gb}); + auto engine_ref = Engine(Operator::delcgtga_delcgtgb, max_nprim, max_l, 0, + eps, params_t{ga, gb}); + const auto ints = compute_all(engine); + const auto ints_ref = compute_all(engine_ref); + REQUIRE(ints.size() == ints_ref.size()); + for (std::size_t i = 0; i != ints.size(); ++i) + REQUIRE(ints[i] == Approx(ints_ref[i]).margin(1e-12)); + } } // see https://github.com/evaleev/libint/issues/216 diff --git a/include/libint2/engine.h b/include/libint2/engine.h index 1ce0f189d..bbb04274a 100644 --- a/include/libint2/engine.h +++ b/include/libint2/engine.h @@ -225,9 +225,11 @@ enum class Operator { stg_x_coulomb, /// alias for Operator::stg_x_coulomb yukawa = stg_x_coulomb, - /// (Delta . contracted Gaussian geminal a) . (Delta . contracted Gaussian - /// geminal b), i.e. the two-geminal generalization of delcgtg2; params are - /// std::array{a, b} + /// \f$ (\nabla_1 g_a) \cdot (\nabla_1 g_b) \f$ , where \f$ g_a \f$ and + /// \f$ g_b \f$ are contracted Gaussian geminals, i.e. the two-geminal + /// generalization of delcgtg2 (delcgtga_delcgtgb with \f$ g_a = g_b = g \f$ + /// is equivalent to delcgtg2 with \f$ g \f$ ); params (no defaults, must be + /// provided) are std::array{g_a, g_b} delcgtga_delcgtgb, // do not modify this invalid = -1, @@ -459,8 +461,8 @@ template <> struct operator_traits : public detail::cgtg_operator_traits<2> { /// the two geminals a and b + /// N.B. like delcgtg2 there are no default params typedef std::array oper_params_type; - static oper_params_type default_params() { return oper_params_type{}; } }; template <> @@ -803,6 +805,11 @@ class Engine { oper_ = new_oper; braket_ = default_braket(oper_); params_ = default_params(oper_); + // discard the core params of the previous operator; N.B. not using + // init_core_ints_params since it throws for operators without default + // params (e.g. delcgtg2), computing their integrals will throw until + // set_params() is called + core_ints_params_ = params_; initialize(); core_eval_pack_ = make_core_eval_pack(oper_); // must follow initialize() // to ensure default diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index c94c92838..314480617 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -70,31 +70,31 @@ typename std::remove_all_extents::type* to_ptr1(T (&a)[N]) { /// These MUST appear in the same order as in Operator. /// You must also update BOOST_PP_NBODY_OPERATOR_LAST_ONEBODY_INDEX when you add /// one-body ints -#define BOOST_PP_NBODY_OPERATOR_LIST \ - (overlap, /* overlap */ \ - (kinetic, /* kinetic */ \ - (elecpot, /* nuclear */ \ - (elecpot, /* erf_nuclear */ \ - (elecpot, /* erfc_nuclear */ \ - (elecpot, /* erfx_nuclear */ \ - (elecpot, /* q_gau */ \ - (1emultipole, /* emultipole1 */ \ - (2emultipole, /* emultipole2 */ \ - (3emultipole, /* emultipole3 */ \ - (sphemultipole, /* sphemultipole */ \ - (opVop, /* opVop */ \ - (opVop, /* op_q_gau_op */ \ - (eri, /* delta */ \ - (eri, /* coulomb */ \ - (eri, /* cgtg */ \ - (eri, /* cgtg_x_coulomb */ \ - (eri, /* delcgtg2 */ \ - (eri, /* r12 */ \ - (eri, /* erf_coulomb */ \ - (eri, /* erfc_coulomb */ \ - (eri, /* erfx_coulomb */ \ - (eri, /* stg */ \ - (eri, /* yukawa */ \ +#define BOOST_PP_NBODY_OPERATOR_LIST \ + (overlap, /* overlap */ \ + (kinetic, /* kinetic */ \ + (elecpot, /* nuclear */ \ + (elecpot, /* erf_nuclear */ \ + (elecpot, /* erfc_nuclear */ \ + (elecpot, /* erfx_nuclear */ \ + (elecpot, /* q_gau */ \ + (1emultipole, /* emultipole1 */ \ + (2emultipole, /* emultipole2 */ \ + (3emultipole, /* emultipole3 */ \ + (sphemultipole, /* sphemultipole */ \ + (opVop, /* opVop */ \ + (opVop, /* op_q_gau_op */ \ + (eri, /* delta */ \ + (eri, /* coulomb */ \ + (eri, /* cgtg */ \ + (eri, /* cgtg_x_coulomb */ \ + (eri, /* delcgtg2 */ \ + (eri, /* r12 */ \ + (eri, /* erf_coulomb */ \ + (eri, /* erfc_coulomb */ \ + (eri, /* erfx_coulomb */ \ + (eri, /* stg */ \ + (eri, /* yukawa */ \ (eri, /* delcgtga_delcgtgb */ \ BOOST_PP_NIL))))))))))))))))))))))))) @@ -908,14 +908,30 @@ __libint2_engine_inline any Engine::make_core_eval_pack(Operator oper) const { } __libint2_engine_inline void Engine::init_core_ints_params(const any& params) { + // geminal exponents must be nonnegative, else the integrals diverge + auto validate_geminal = [](const ContractedGaussianGeminal& g, + const char* oper_name) { + for (const auto& p : g) + if (p.first < 0) + throw std::logic_error( + "Engine::set_params(Operator::" + std::string(oper_name) + + "): geminal exponent " + std::to_string(p.first) + " is negative"); + }; + if (oper_ == Operator::delcgtg2) { // [g12,[- \Del^2, g12] = 2 (\Del g12) \cdot (\Del g12) // (\Del exp(-a r_12^2) \cdot (\Del exp(-b r_12^2) = 4 a b (r_{12}^2 exp(- // (a+b) r_{12}^2) ) // i.e. need to scale each coefficient by 4 a b - const auto& oparams = - any_cast::oper_params_type&>( - params); + const auto* oparams_ptr = + any_cast::oper_params_type>( + ¶ms); + if (oparams_ptr == nullptr) + throw std::logic_error( + "Engine::set_params(Operator::delcgtg2): params must be provided as " + "a ContractedGaussianGeminal"); + const auto& oparams = *oparams_ptr; + validate_geminal(oparams, "delcgtg2"); const auto ng = oparams.size(); operator_traits::oper_params_type core_ints_params; core_ints_params.reserve(ng * (ng + 1) / 2); @@ -948,10 +964,17 @@ __libint2_engine_inline void Engine::init_core_ints_params(const any& params) { // g_b = \sum_j d_j exp(-b_j r_12^2): // = \sum_ij 4 a_i b_j c_i d_j r_{12}^2 exp(-(a_i+b_j) r_{12}^2) // i.e. product of the two geminals, each coefficient scaled by 4 a_i b_j - const auto& oparams = any_cast::oper_params_type&>(params); - const auto& ga = oparams[0]; - const auto& gb = oparams[1]; + const auto* oparams_ptr = any_cast< + operator_traits::oper_params_type>( + ¶ms); + if (oparams_ptr == nullptr) + throw std::logic_error( + "Engine::set_params(Operator::delcgtga_delcgtgb): params must be " + "provided as a std::array"); + const auto& ga = (*oparams_ptr)[0]; + const auto& gb = (*oparams_ptr)[1]; + validate_geminal(ga, "delcgtga_delcgtgb"); + validate_geminal(gb, "delcgtga_delcgtgb"); operator_traits::oper_params_type core_ints_params; core_ints_params.reserve(ga.size() * gb.size()); for (const auto& pa : ga) @@ -1518,23 +1541,21 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( core_ints_params_); core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params); } break; - case Operator::delcgtg2: { - const auto& core_eval_ptr = - any_cast&>(core_eval_pack_) - .first(); - const auto& core_ints_params = - any_cast::oper_params_type&>( - core_ints_params_); - core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params); - } break; + // N.B. delcgtga_delcgtgb has the same core evaluator as + // delcgtg2, and for both core_ints_params_ holds the product + // geminal, see init_core_ints_params + case Operator::delcgtg2: case Operator::delcgtga_delcgtgb: { - // core_ints_params_ holds the product geminal, see - // init_core_ints_params + static_assert( + std::is_same< + detail::core_eval_pack_type, + detail::core_eval_pack_type< + Operator::delcgtga_delcgtgb>>::value, + "delcgtg2 and delcgtga_delcgtgb must share the core " + "evaluator"); const auto& core_eval_ptr = any_cast&>(core_eval_pack_) + Operator::delcgtg2>&>(core_eval_pack_) .first(); const auto& core_ints_params = any_cast Date: Sun, 25 Oct 2026 11:17:28 -0400 Subject: [PATCH 4/4] Engine: instantiate enforce_params_type for geminal params in the compiled library With LIBINT2_DOES_NOT_INLINE_ENGINE, enforce_params_type was only instantiated for empty params, so constructing an Engine for delcgtg2 or delcgtga_delcgtgb (which have no default params and must be constructed with params) failed to link against libint2-cxx. --- include/libint2/engine.impl.h | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index 314480617..95a54b16b 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -2273,6 +2273,18 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( template any Engine::enforce_params_type( Operator oper, const Engine::empty_pod& params, bool throw_if_wrong_type); +// operators without default params (e.g. delcgtg2) can only be constructed +// with params, hence instantiate for their params types +template any Engine::enforce_params_type( + Operator oper, const ContractedGaussianGeminal& params, + bool throw_if_wrong_type); +template any Engine::enforce_params_type< + operator_traits::oper_params_type>( + Operator oper, + const operator_traits::oper_params_type& + params, + bool throw_if_wrong_type); + template const Engine::target_ptr_vec& Engine::compute( const Shell& first_shell, const Shell&);