Relativistic 2e integrals

cmake/modules/int_am.cmake

@@ -262,7 +262,7 @@ macro(process_integrals_class class)
             if (LIBINT2_${class}_OPT_AM EQUAL -1)
                 set(LIBINT_${class}_OPT_AM "")
             else()
-                set($LIBINT_{class}_OPT_AM ${LIBINT2_${class}_OPT_AM})
+                set(LIBINT_${class}_OPT_AM ${LIBINT2_${class}_OPT_AM})
             endif()
         endif()
         if (LIBINT_OPT_AM_LIST)
@@ -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.h

@@ -22,18 +22,19 @@
 #define _libint2_header_
 
 #define LIBINT_T_SS_EREP_SS(mValue) \
-  _aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_##mValue
+  _aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_##mValue
 #define LIBINT_T_SS_Km1G12_SS(mValue) \
-  _aB_s___0__s___1___r12_minus_1_g12_s___0__s___1___Ab__up_##mValue
+  _aB_s____0__s____1___r12_minus_1_g12_s____0__s____1___Ab__up_##mValue
 #define LIBINT_T_SS_K0G12_SS_0 \
-  _aB_s___0__s___1___r12_0_g12_s___0__s___1___Ab__up_0
+  _aB_s____0__s____1___r12_0_g12_s____0__s____1___Ab__up_0
 #define LIBINT_T_SS_K2G12_SS_0 \
-  _aB_s___0__s___1___r12_2_g12_s___0__s___1___Ab__up_0
+  _aB_s____0__s____1___r12_2_g12_s____0__s____1___Ab__up_0
 #define LIBINT_T_SS_K4G12_SS_0 \
-  _aB_s___0__s___1___r12_4_g12_s___0__s___1___Ab__up_0
-#define LIBINT_T_S_OVERLAP_S _aB_s___0___Overlap_s___0___Ab__up_
-#define LIBINT_T_S_KINETIC_S _aB_s___0___Kinetic_s___0___Ab__up_
-#define LIBINT_T_S_ELECPOT_S(mValue) _aB_s___0___ElecPot_s___0___Ab__up_##mValue
+  _aB_s____0__s____1___r12_4_g12_s____0__s____1___Ab__up_0
+#define LIBINT_T_S_OVERLAP_S _aB_s____0___Overlap_s____0___Ab__up_
+#define LIBINT_T_S_KINETIC_S _aB_s____0___Kinetic_s____0___Ab__up_
+#define LIBINT_T_S_ELECPOT_S(mValue) \
+  _aB_s____0___ElecPot_s____0___Ab__up_##mValue
 
 #include <libint2/util/configuration.h>
 #include <libint2/util/generated/libint2_params.h>

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_RKB_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_RKB_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,13 @@ 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,
+  /// (2-body) \f$ (σ.p_{b1})(σ.p_{b2}) r_{12}^{-1} (σ.p_{k1})(σ.p_{k2})\f$
+  /// where b1 & b2 are centers of bra1 and bra2 and k1  & k2 are centers of
+  /// ket1 and ket2, respectively
+  opop_coulomb_opop,
   /// contracted Gaussian geminal
   cgtg,
   /// contracted Gaussian geminal times Coulomb
@@ -246,6 +253,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 +354,22 @@ 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;
+};
+template <>
+struct operator_traits<Operator::opop_coulomb_opop>
+    : public operator_traits<Operator::coulomb> {
+  /// 16 components: tensor product of two independent spin-space quaternions
+  /// index = 4 * bra_spin + ket_spin, where spin in {S=0, X=1, Y=2, Z=3}
+  static constexpr auto nopers = 16;
+  static constexpr auto intrinsic_deriv_order = 4;
+};
+
 namespace detail {
 template <int K>
 struct cgtg_operator_traits : public detail::default_operator_traits {
@@ -839,16 +863,16 @@ class Engine {
       const Shell& ket2, const ShellPair* spbra, const ShellPair* spket);
 
   // clang-format off
-  /** this specifies target precision for computing the integrals, i.e.
-   *  the target absolute (i.e., not relative) error of the integrals.
-   *  It is used to screen out primitive integrals. For some screening
-   *  methods precision can be almost guaranteed (due to finite precision
-   *  of the precomputed interpolation tables used to evaluate the core integrals
-   *  it is not in general possible to guarantee precision rigorously).
-   *
-   *  @param[in] prec the target precision
-   *  @sa ScreeningMethod
-   */
+ /** this specifies target precision for computing the integrals, i.e.
+  *  the target absolute (i.e., not relative) error of the integrals.
+  *  It is used to screen out primitive integrals. For some screening
+  *  methods precision can be almost guaranteed (due to finite precision
+  *  of the precomputed interpolation tables used to evaluate the core integrals
+  *  it is not in general possible to guarantee precision rigorously).
+  *
+  *  @param[in] prec the target precision
+  *  @sa ScreeningMethod
+  */
   // clang-format on
   Engine& set_precision(scalar_type prec) {
     if (prec <= 0.) {