diff --git a/export/tests/unit/test-2body.cc b/export/tests/unit/test-2body.cc index fd602a910..d25bcc243 100644 --- a/export/tests/unit/test-2body.cc +++ b/export/tests/unit/test-2body.cc @@ -122,6 +122,141 @@ 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) { + 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(), target, target + setsize); + } + } + } + } + } + return ints; + }; + + 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 + 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(params_t{gb, ga}); + const auto ints_ba = compute_all(engine); + 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); + } + } + + 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 TEST_CASE_METHOD(libint2::unit::DefaultFixture, "bra-ket permutation", "[engine][2-body]") { diff --git a/include/libint2/engine.h b/include/libint2/engine.h index f07e8096c..bbb04274a 100644 --- a/include/libint2/engine.h +++ b/include/libint2/engine.h @@ -225,6 +225,12 @@ enum class Operator { stg_x_coulomb, /// alias for Operator::stg_x_coulomb yukawa = stg_x_coulomb, + /// \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, // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!keep this @@ -232,7 +238,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 +457,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 + /// N.B. like delcgtg2 there are no default params + typedef std::array oper_params_type; +}; template <> struct operator_traits @@ -792,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 f2f7cedcd..95a54b16b 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -70,32 +70,33 @@ 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 */ \ - BOOST_PP_NIL)))))))))))))))))))))))) +#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))))))))))))))))))))))))) #define BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE \ BOOST_PP_MAKE_TUPLE(BOOST_PP_LIST_SIZE(BOOST_PP_NBODY_OPERATOR_LIST)) @@ -907,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); @@ -942,6 +959,32 @@ __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_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) + 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; } @@ -1498,7 +1541,18 @@ __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: { + // 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: { + 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_) @@ -2219,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&);