Add operator kernel delcgtga_delcgtgb

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<Shell> 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<std::pair<double, double>>;
+  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<double>::epsilon();
+  using params_t = std::array<geminal_t, 2>;
+
+  // computes all shell quartets of obs (all targets, i.e. all derivatives)
+  // with engine e
+  auto compute_all = [&](Engine &e) {
+    std::vector<double> 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]") {

include/libint2/engine.h

@@ -225,14 +225,20 @@ 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<ContractedGaussianGeminal,2>{g_a, g_b}
+  delcgtga_delcgtgb,
   // do not modify this
   invalid = -1,
   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!keep this
   // updated!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   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<Operator::cgtg_x_coulomb>
 template <>
 struct operator_traits<Operator::delcgtg2>
     : public detail::cgtg_operator_traits<2> {};
+template <>
+struct operator_traits<Operator::delcgtga_delcgtgb>
+    : public detail::cgtg_operator_traits<2> {
+  /// the two geminals a and b
+  /// N.B. like delcgtg2 there are no default params
+  typedef std::array<ContractedGaussianGeminal, 2> oper_params_type;
+};
 
 template <>
 struct operator_traits<Operator::delta>
@@ -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

include/libint2/engine.impl.h

@@ -70,32 +70,33 @@ typename std::remove_all_extents<T>::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<const operator_traits<Operator::delcgtg2>::oper_params_type&>(
-            params);
+    const auto* oparams_ptr =
+        any_cast<operator_traits<Operator::delcgtg2>::oper_params_type>(
+            &params);
+    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<Operator::delcgtg2>::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<Operator::delcgtga_delcgtgb>::oper_params_type>(
+        &params);
+    if (oparams_ptr == nullptr)
+      throw std::logic_error(
+          "Engine::set_params(Operator::delcgtga_delcgtgb): params must be "
+          "provided as a std::array<ContractedGaussianGeminal,2>");
+    const auto& ga = (*oparams_ptr)[0];
+    const auto& gb = (*oparams_ptr)[1];
+    validate_geminal(ga, "delcgtga_delcgtgb");
+    validate_geminal(gb, "delcgtga_delcgtgb");
+    operator_traits<Operator::cgtg>::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<Operator::delcgtg2>,
+                          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<const detail::core_eval_pack_type<
                           Operator::delcgtg2>&>(core_eval_pack_)
@@ -2219,6 +2273,18 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
 template any Engine::enforce_params_type<Engine::empty_pod>(
     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<ContractedGaussianGeminal>(
+    Operator oper, const ContractedGaussianGeminal& params,
+    bool throw_if_wrong_type);
+template any Engine::enforce_params_type<
+    operator_traits<Operator::delcgtga_delcgtgb>::oper_params_type>(
+    Operator oper,
+    const operator_traits<Operator::delcgtga_delcgtgb>::oper_params_type&
+        params,
+    bool throw_if_wrong_type);
+
 template const Engine::target_ptr_vec& Engine::compute<Shell>(
     const Shell& first_shell, const Shell&);