evaleev · kshitij-05 · Feb 5, 2026 · Feb 7, 2026 · Feb 9, 2026 · Feb 10, 2026
diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml
@@ -63,6 +63,8 @@ jobs:
         -D LIBINT2_ERI3_MAX_AM=3;2
         -D LIBINT2_ENABLE_ERI=1
         -D LIBINT2_ENABLE_ERI3=1
+        -D LIBINT2_ENABLE_RKB_ERI3=0
+        -D LIBINT2_RKB_ERI3_MAX_AM=2
         -D LIBINT2_ENABLE_ONEBODY=1
         -D LIBINT2_DISABLE_ONEBODY_PROPERTY_DERIVS=ON
         -D LIBINT2_MULTIPOLE_MAX_ORDER=2

diff --git a/CMakeLists.txt b/CMakeLists.txt
diff --git a/INSTALL.md b/INSTALL.md
@@ -169,6 +169,8 @@ Use combined targets like `cmake --target check install` to avoid some unnecessa
 * `LIBINT2_ENABLE_ERI` - G - Compile with support for up to N-th derivatives of 4-center electron repulsion integrals. Use -1 for OFF. [Default=0]
 * `LIBINT2_ENABLE_ERI3` - G - Compile with support for up to N-th derivatives of 3-center electron repulsion integrals. Use -1 for OFF. [Default=-1]
 * `LIBINT2_ENABLE_ERI2` - G - Compile with support for up to N-th derivatives of 2-center electron repulsion integrals. Use -1 for OFF. [Default=-1]
+* `LIBINT2_ENABLE_RKB_ERI` - G - Compile with support for up to N-th derivatives of relativistic restricted kinetic balance (RKB) 4-center electron repulsion integrals. Use -1 for OFF. [Default=0]
+* `LIBINT2_ENABLE_RKB_ERI3` - G - Compile with support for up to N-th derivatives of relativistic restricted kinetic balance (RKB) 3-center electron repulsion integrals. Use -1 for OFF. σ·p acts on the two paired (AO) centers; the unpaired/fitting (density-fitting/auxiliary) center is a spectator. [Default=-1]
 * `LIBINT2_ENABLE_G12` - G - Compile with support for N-th derivatives of MP2-F12 energies with Gaussian factors. Use -1 for OFF. [Default=-1]
 * `LIBINT2_ENABLE_G12DKH` - G - Compile with support for N-th derivatives of DKH-MP2-F12 energies with Gaussian factors. Use -1 for OFF. [Default=-1]
 
@@ -255,6 +257,12 @@ Note that options, docs, and CMake components are focused on the C++ interface,
 * `LIBINT2_ERI2_OPT_AM` - G - Optimize 2-center ERIs maximally for up to angular momentum N (N <= max-am). Can specify values for each derivative level as a semicolon-separated string. [Default=-1 -> `LIBINT2_OPT_AM`]
 * `LIBINT2_ERI2_PURE_SH` - G - Assume the 2-center ERIs will be transformed to pure solid harmonics. [Default=OFF]
 
+* `LIBINT2_RKB_ERI_MAX_AM` - G - Support RKB 4-center ERIs for Gaussians of angular momentum up to N. Can specify values for each derivative level as a semicolon-separated string. [Default=-1 -> `LIBINT2_MAX_AM`]
+* `LIBINT2_RKB_ERI_OPT_AM` - G - Optimize RKB 4-center ERIs maximally for up to angular momentum N (N <= max-am). Can specify values for each derivative level as a semicolon-separated string. [Default=-1 -> `LIBINT2_OPT_AM`]
+
+* `LIBINT2_RKB_ERI3_MAX_AM` - G - Support RKB 3-center ERIs for Gaussians of angular momentum up to N. Can specify values for each derivative level as a semicolon-separated string. This option sets the AM of the density-fitting (auxiliary) basis on the single fitting center; the paired (AO) centers, on which σ·p acts, use LIBINT2_MAX_AM. [Default=-1 -> `LIBINT2_MAX_AM`]
+* `LIBINT2_RKB_ERI3_OPT_AM` - G - Optimize RKB 3-center ERIs maximally for up to angular momentum N (N <= max-am). Can specify values for each derivative level as a semicolon-separated string. [Default=-1 -> `LIBINT2_OPT_AM`]
+
 * `LIBINT2_G12_MAX_AM` - G - Support integrals for G12 methods of angular momentum up to N. No specification with per-derivative list. [Default=-1 -> `LIBINT2_MAX_AM`]
 * `LIBINT2_G12_OPT_AM` - G - Optimize G12 integrals for up to angular momentum N (N <= max-am). No specification with per-derivative list. [Default=-1 `LIBINT2_OPT_AM`]
 

diff --git a/cmake/modules/int_am.cmake b/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,8 @@ endmacro()
 
 process_integrals_class(ONEBODY)
 process_integrals_class(ERI)
+process_integrals_class(RKB_ERI)
+process_integrals_class(RKB_ERI3)
 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
@@ -380,6 +382,7 @@ add_feature_info(
 set(Libint2_ERI_COMPONENTS "")
 set(_amlist "")
 set(_eri3_impure_sh "")
+set(_rkb_eri3_impure_sh "")
 set(_eri2_impure_sh "")
 add_feature_info(
   "integral class MULTIPOLE derivative 0"
@@ -396,7 +399,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;RKB_ERI3;ERI3;ERI2;G12;G12DKH)
     if((_cls STREQUAL G12) OR (_cls STREQUAL G12DKH))
         add_feature_info(
           "integral class ${_cls}"
@@ -433,6 +436,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             "rkberi_${_am${_l}}${_am${_l}}${_am${_l}}${_am${_l}}_d${_d}")
                 endif()
             endforeach()
             if (_cls STREQUAL "ERI3")
@@ -445,6 +450,19 @@ foreach(_cls ONEBODY;ERI;ERI3;ERI2;G12;G12DKH)
                     endforeach()
                 endforeach()
             endif()
+            if (_cls STREQUAL "RKB_ERI3")
+                # Mirror ERI3 component naming: fitting (single) center on bra,
+                # paired (AO) centers on ket where σ·p acts. Paired-center AM
+                # tracks LIBINT_MAX_AM via _eri3_candidate0_d${_d}.
+                foreach(_lfit RANGE ${LIBINT_HARD_MIN_AM} ${_candidate_${_cls}_d${_d}})  # LIBINT_RKB_ERI3_MAX_AM[_LIST], fitting
+                    foreach(_lpr RANGE ${LIBINT_HARD_MIN_AM} ${_eri3_candidate0_d${_d}})  # LIBINT_MAX_AM[_LIST], paired
+                        if (_lfit GREATER_EQUAL _lpr)
+                            list(APPEND _amlist     "rkberi_${_am${_lpr}}${_am${_lpr}}${_AM${_lfit}}_d${_d}")
+                            list(APPEND _pureamlist "rkberi_${_am${_lpr}}${_am${_lpr}}${_am${_lfit}}_d${_d}")
+                        endif()
+                    endforeach()
+                endforeach()
+            endif()
             list(REVERSE _amlist)
             list(APPEND Libint2_ERI_COMPONENTS "${_amlist}")
             message(VERBOSE "setting components ${_amlist}")
@@ -453,11 +471,16 @@ foreach(_cls ONEBODY;ERI;ERI3;ERI2;G12;G12DKH)
                 list(APPEND _eri2_impure_sh "${_pureamlist}")
             elseif (_cls STREQUAL "ERI3")
                 list(APPEND _eri3_impure_sh "${_pureamlist}")
+            elseif (_cls STREQUAL "RKB_ERI3")
+                list(APPEND _rkb_eri3_impure_sh "${_pureamlist}")
             endif()
         endforeach()
         if ((_cls STREQUAL "ERI3") AND NOT LIBINT2_ERI3_PURE_SH)
             list(APPEND Libint2_ERI_COMPONENTS "${_eri3_impure_sh}")
             message(VERBOSE "setting components ${_eri3_impure_sh}")
+        elseif ((_cls STREQUAL "RKB_ERI3") AND NOT LIBINT2_RKB_ERI3_PURE_SH)
+            list(APPEND Libint2_ERI_COMPONENTS "${_rkb_eri3_impure_sh}")
+            message(VERBOSE "setting components ${_rkb_eri3_impure_sh}")
         elseif ((_cls STREQUAL "ERI2") AND NOT LIBINT2_ERI2_PURE_SH)
             list(APPEND Libint2_ERI_COMPONENTS "${_eri2_impure_sh}")
             message(VERBOSE "setting components ${_eri2_impure_sh}")

diff --git a/export/fortran/libint_f.F90 b/export/fortran/libint_f.F90
@@ -520,88 +520,88 @@ SUBROUTINE compute_eri_f(contrdepth, deriv_order, am1, c1, alpha1, A, &
 
                   am_tot = am + deriv_order
 
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_0
-                  IF (am_tot >= 0) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_0(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_0
+                  IF (am_tot >= 0) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_0(1) = &
                      C1234*pfac*F(1, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_1
-                  IF (am_tot >= 1) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_1(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_1
+                  IF (am_tot >= 1) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_1(1) = &
                      C1234*pfac*F(2, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_2
-                  IF (am_tot >= 2) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_2(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_2
+                  IF (am_tot >= 2) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_2(1) = &
                      C1234*pfac*F(3, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_3
-                  IF (am_tot >= 3) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_3(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_3
+                  IF (am_tot >= 3) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_3(1) = &
                      C1234*pfac*F(4, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_4
-                  IF (am_tot >= 4) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_4(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_4
+                  IF (am_tot >= 4) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_4(1) = &
                      C1234*pfac*F(5, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_5
-                  IF (am_tot >= 5) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_5(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_5
+                  IF (am_tot >= 5) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_5(1) = &
                      C1234*pfac*F(6, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_6
-                  IF (am_tot >= 6) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_6(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_6
+                  IF (am_tot >= 6) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_6(1) = &
                      C1234*pfac*F(7, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_7
-                  IF (am_tot >= 7) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_7(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_7
+                  IF (am_tot >= 7) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_7(1) = &
                      C1234*pfac*F(8, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_8
-                  IF (am_tot >= 8) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_8(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_8
+                  IF (am_tot >= 8) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_8(1) = &
                      C1234*pfac*F(9, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_9
-                  IF (am_tot >= 9) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_9(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_9
+                  IF (am_tot >= 9) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_9(1) = &
                      C1234*pfac*F(10, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_10
-                  IF (am_tot >= 10) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_10(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_10
+                  IF (am_tot >= 10) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_10(1) = &
                      C1234*pfac*F(11, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_11
-                  IF (am_tot >= 11) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_11(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_11
+                  IF (am_tot >= 11) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_11(1) = &
                      C1234*pfac*F(12, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_12
-                  IF (am_tot >= 12) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_12(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_12
+                  IF (am_tot >= 12) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_12(1) = &
                      C1234*pfac*F(13, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_13
-                  IF (am_tot >= 13) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_13(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_13
+                  IF (am_tot >= 13) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_13(1) = &
                      C1234*pfac*F(14, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_14
-                  IF (am_tot >= 14) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_14(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_14
+                  IF (am_tot >= 14) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_14(1) = &
                      C1234*pfac*F(15, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_15
-                  IF (am_tot >= 15) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_15(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_15
+                  IF (am_tot >= 15) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_15(1) = &
                      C1234*pfac*F(16, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_16
-                  IF (am_tot >= 16) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_16(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_16
+                  IF (am_tot >= 16) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_16(1) = &
                      C1234*pfac*F(17, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_17
-                  IF (am_tot >= 17) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_17(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_17
+                  IF (am_tot >= 17) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_17(1) = &
                      C1234*pfac*F(18, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_18
-                  IF (am_tot >= 18) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_18(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_18
+                  IF (am_tot >= 18) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_18(1) = &
                      C1234*pfac*F(19, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_19
-                  IF (am_tot >= 19) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_19(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_19
+                  IF (am_tot >= 19) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_19(1) = &
                      C1234*pfac*F(20, p1234)
 #endif
-#ifdef LIBINT2_DEFINED__aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_20
-                  IF (am_tot >= 20) erieval(p1234)%f_aB_s___0__s___1___TwoPRep_s___0__s___1___Ab__up_20(1) = &
+#ifdef LIBINT2_DEFINED__aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_20
+                  IF (am_tot >= 20) erieval(p1234)%f_aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_20(1) = &
                      C1234*pfac*F(21, p1234)
 #endif
                ENDDO