Relativistic 2e 3- and 4-center integrals

cmake/modules/int_am.cmake

@@ -357,6 +357,7 @@ endmacro()
 
 process_integrals_class(ONEBODY)
 process_integrals_class(ERI)
+process_integrals_class(RKB_ERI)
 process_integrals_class(ERI3)
 process_integrals_class(ERI2)
 # unlike above, these classes (1) don't do AM_LIST and (2) require value in config.h if enabled
@@ -396,7 +397,7 @@ list(REVERSE _amlist)
 list(APPEND Libint2_ERI_COMPONENTS "${_amlist}")
 message(VERBOSE "setting components ${_amlist}")
 
-foreach(_cls ONEBODY;ERI;ERI3;ERI2;G12;G12DKH)
+foreach(_cls ONEBODY;ERI;RKB_ERI;ERI3;ERI2;G12;G12DKH)
     if((_cls STREQUAL G12) OR (_cls STREQUAL G12DKH))
         add_feature_info(
           "integral class ${_cls}"
@@ -433,6 +434,8 @@ foreach(_cls ONEBODY;ERI;ERI3;ERI2;G12;G12DKH)
                     list(APPEND _amlist         "onebody_${_am${_l}}${_am${_l}}_d${_d}")
                 elseif (_cls STREQUAL "G12")
                     list(APPEND _amlist             "g12_${_am${_l}}${_am${_l}}${_am${_l}}${_am${_l}}_d${_d}")
+                elseif (_cls STREQUAL "RKB_ERI")
+                    list(APPEND _amlist             "rkb_eri_${_am${_l}}${_am${_l}}${_am${_l}}${_am${_l}}_d${_d}")
                 endif()
             endforeach()
             if (_cls STREQUAL "ERI3")

include/libint2/config.h.cmake.in

@@ -71,6 +71,13 @@
 #undef LIBINT_INCLUDE_ERI
 #endif
 
+/* Support ERI derivatives up to this order */
+#define LIBINT_INCLUDE_RKB_ERI @LIBINT_INCLUDE_RKB_ERI@
+#if @LIBINT_INCLUDE_RKB_ERI@ == -1
+#undef LIBINT_INCLUDE_RKB_ERI
+#endif
+
+
 /* Support 3-center ERI derivatives up to this order */
 #define LIBINT_INCLUDE_ERI3 @LIBINT_INCLUDE_ERI3@
 #if @LIBINT_INCLUDE_ERI3@ == -1
@@ -122,6 +129,18 @@
 /* Max optimized AM for ERI and its derivatives */
 #cmakedefine LIBINT_ERI_OPT_AM_LIST "@LIBINT_ERI_OPT_AM_LIST@"
 
+/* Max AM for RKB_ERI (same for all derivatives; if not defined see LIBINT_ERI_MAX_AM_LIST) */
+#cmakedefine LIBINT_RKB_ERI_MAX_AM @LIBINT_RKB_ERI_MAX_AM@
+
+/* Max AM for RKB_ERI and its derivatives */
+#cmakedefine LIBINT_RKB_ERI_MAX_AM_LIST "@LIBINT_RKB_ERI_MAX_AM_LIST@"
+
+/* Max optimized AM for ERI (same for all derivatives; if not defined see LIBINT_ERI_OPT_AM_LIST) */
+#cmakedefine LIBINT_RKB_ERI_OPT_AM @LIBINT_RKB_ERI_OPT_AM@
+
+/* Max optimized AM for ERI and its derivatives */
+#cmakedefine LIBINT_RKB_ERI_OPT_AM_LIST "@LIBINT_RKB_ERI_OPT_AM_LIST@"
+
 /* Max AM for 3-center ERI (same for all derivatives; if not defined see LIBINT_ERI3_MAX_AM_LIST) */
 #cmakedefine LIBINT_ERI3_MAX_AM @LIBINT_ERI3_MAX_AM@
 

include/libint2/cxxapi.h

@@ -35,9 +35,9 @@
 
 #if !defined(LIBINT_INCLUDE_ONEBODY) ||                              \
     !(defined(LIBINT_INCLUDE_ERI) || defined(LIBINT_INCLUDE_ERI3) || \
-      defined(LIBINT_INCLUDE_ERI2))
+      defined(LIBINT_INCLUDE_ERI2) || defined(LIBINT_INCLUDE_RKB_ERI))
 #error \
-    "C++ API is only supported if both 1-body and some (eri, eri3, eri2) 2-body integrals are enabled"
+    "C++ API is only supported if both 1-body and some (eri, eri3, eri2, rkb_eri) 2-body integrals are enabled"
 #endif
 
 #include <libint2/atom.h>

include/libint2/engine.h

@@ -153,6 +153,9 @@ enum class Operator {
   coulomb,
   /// alias for Operator::coulomb
   r12_m1 = coulomb,
+  /// (2-body) \f$ r_{12}^{-1} (σ.p_{k1})(σ.p_{k2})\f$ where k1  & k2 are
+  /// centers of ket1 and ket2, respectively
+  coulomb_opop,
   /// contracted Gaussian geminal
   cgtg,
   /// contracted Gaussian geminal times Coulomb
@@ -246,6 +249,7 @@ struct operator_traits<Operator::nuclear>
   typedef const libint2::FmEval_Reference<scalar_type> core_eval_type;
 #endif
 };
+
 template <>
 struct operator_traits<Operator::opVop>
     : public operator_traits<Operator::nuclear> {
@@ -346,6 +350,14 @@ struct operator_traits<Operator::coulomb>
   typedef const libint2::FmEval_Reference<scalar_type> core_eval_type;
 #endif
 };
+
+template <>
+struct operator_traits<Operator::coulomb_opop>
+    : public operator_traits<Operator::coulomb> {
+  static constexpr auto nopers = 4;
+  static constexpr auto intrinsic_deriv_order = 2;
+};
+
 namespace detail {
 template <int K>
 struct cgtg_operator_traits : public detail::default_operator_traits {

include/libint2/engine.impl.h

@@ -70,30 +70,31 @@ 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 */   \
-        (1emultipole,        /* emultipole1 */    \
-         (2emultipole,       /* emultipole2 */    \
-          (3emultipole,      /* emultipole3 */    \
-           (sphemultipole,   /* sphemultipole */  \
-            (opVop,          /* opVop */          \
-             (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 */   \
+        (1emultipole,         /* emultipole1 */    \
+         (2emultipole,        /* emultipole2 */    \
+          (3emultipole,       /* emultipole3 */    \
+           (sphemultipole,    /* sphemultipole */  \
+            (opVop,           /* opVop */          \
+             (eri,            /* delta */          \
+              (eri,           /* coulomb */        \
+               (eri,          /* coulomb_opop */   \
+                (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_INDEX_TUPLE \
   BOOST_PP_MAKE_TUPLE(BOOST_PP_LIST_SIZE(BOOST_PP_NBODY_OPERATOR_LIST))
@@ -663,23 +664,23 @@ __libint2_engine_inline void Engine::initialize(size_t max_nprim) {
   // validate braket
 #ifndef LIBINT_INCLUDE_ONEBODY
   assert(braket_ != BraKet::x_x &&
-         "this braket type not supported by the library; give --enable-1body "
-         "to configure");
+         "this braket type not supported by the library; configure with "
+         "-DLIBINT_INCLUDE_ONEBODY >= 0");
 #endif
 #ifndef LIBINT_INCLUDE_ERI
   assert(braket_ != BraKet::xx_xx &&
-         "this braket type not supported by the library; give --enable-eri to "
-         "configure");
+         "this braket type not supported by the library; configure with "
+         "-DLIBINT_INCLUDE_ERI >= 0");
 #endif
 #ifndef LIBINT_INCLUDE_ERI3
   assert((braket_ != BraKet::xs_xx && braket_ != BraKet::xx_xs) &&
-         "this braket type not supported by the library; give --enable-eri3 to "
-         "configure");
+         "this braket type not supported by the library; configure with "
+         "-DLIBINT_INCLUDE_ERI3 >= 0");
 #endif
 #ifndef LIBINT_INCLUDE_ERI2
   assert(braket_ != BraKet::xs_xs &&
-         "this braket type not supported by the library; give --enable-eri2 to "
-         "configure");
+         "this braket type not supported by the library; configure with "
+         "-DLIBINT_INCLUDE_ERI2 >= 0");
 #endif
 
   // make sure it's no default initialized
@@ -1421,7 +1422,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_;
+            const auto mmax = l + deriv_order_ + intrinsic_deriv_order();
 
             if (!skip_core_ints) {
               switch (oper_) {
@@ -1432,6 +1433,13 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
                           .first();
                   core_eval_ptr->eval(gm_ptr, T, mmax);
                 } break;
+                case Operator::coulomb_opop: {
+                  const auto& core_eval_ptr =
+                      any_cast<const detail::core_eval_pack_type<
+                          Operator::coulomb_opop>&>(core_eval_pack_)
+                          .first();
+                  core_eval_ptr->eval(gm_ptr, T, mmax);
+                } break;
                 case Operator::cgtg_x_coulomb: {
                   const auto& core_eval_ptr =
                       any_cast<const detail::core_eval_pack_type<