Runtime σ.p integral evaluation via SigmaPEngine (include/libint2/sigma_p.h)

Author: kshitij-05

export/tests/unit/test-2body.cc

@@ -742,6 +742,109 @@ TEST_CASE("RKB Coulomb integrals", "[engine][2-body]") {
     }
   }
 
+  SECTION("Gaunt-type (SL|SL)") {
+    // Need standard ERI support at max_l+1 for AM-shifted integrals
+    if (LIBINT2_MAX_AM_eri < max_l + 1) return;
+
+    libint2::SigmaPEngine sp_engine_gaunt(max_nprim, max_l, BraKet::xx_xx);
+
+    const auto nshell = obs.size();
+    for (int s0 = 0; s0 != nshell; ++s0) {
+      for (int s1 = 0; s1 != nshell; ++s1) {
+        for (int s2 = 0; s2 != nshell; ++s2) {
+          for (int s3 = 0; s3 != nshell; ++s3) {
+            // ([σ.p a] b | [σ.p c] d) — 4 quaternion components
+            const auto &results = sp_engine_gaunt.compute(
+                sigma_p(obs[s0]), obs[s1], sigma_p(obs[s2]), obs[s3]);
+            assert(results.size() == 4);
+
+            LIBINT2_REF_REALTYPE Aref[3], Bref[3], Cref[3], Dref[3];
+            for (int i = 0; i < 3; ++i) Aref[i] = obs[s0].O[i];
+            for (int i = 0; i < 3; ++i) Bref[i] = obs[s1].O[i];
+            for (int i = 0; i < 3; ++i) Cref[i] = obs[s2].O[i];
+            for (int i = 0; i < 3; ++i) Dref[i] = obs[s3].O[i];
+
+            int ijkl = 0;
+
+            int l0, m0, n0;
+            FOR_CART(l0, m0, n0, obs[s0].contr[0].l)
+            int l1, m1, n1;
+            FOR_CART(l1, m1, n1, obs[s1].contr[0].l)
+            int l2, m2, n2;
+            FOR_CART(l2, m2, n2, obs[s2].contr[0].l)
+            int l3, m3, n3;
+            FOR_CART(l3, m3, n3, obs[s3].contr[0].l)
+
+            std::array<LIBINT2_REF_REALTYPE, 4> ref{};
+            ref.fill(0.0);
+
+            for (uint p0 = 0; p0 < obs[s0].nprim(); p0++) {
+              for (uint p1 = 0; p1 < obs[s1].nprim(); p1++) {
+                for (uint p2 = 0; p2 < obs[s2].nprim(); p2++) {
+                  for (uint p3 = 0; p3 < obs[s3].nprim(); p3++) {
+                    const LIBINT2_REF_REALTYPE c0123 =
+                        obs[s0].contr[0].coeff[p0] *
+                        obs[s1].contr[0].coeff[p1] *
+                        obs[s2].contr[0].coeff[p2] * obs[s3].contr[0].coeff[p3];
+
+                    constexpr int X = 0, Y = 1, Z = 2;
+                    // D(r1, r2) = ∂²ERI/∂A_r1 ∂C_r2
+                    auto D = [&](int r1, int r2) {
+                      typedef std::array<unsigned int, 12> didx;
+                      didx d = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+                      d[0 + r1] = 1;  // ∂/∂A_r1 (bra center 1)
+                      d[6 + r2] = 1;  // ∂/∂C_r2 (ket center 1)
+                      return eri(d.data(), l0, m0, n0, obs[s0].alpha[p0], Aref,
+                                 l1, m1, n1, obs[s1].alpha[p1], Bref, l2, m2,
+                                 n2, obs[s2].alpha[p2], Cref, l3, m3, n3,
+                                 obs[s3].alpha[p3], Dref, 0);
+                    };
+
+                    // Quaternion components from (σ.p_a)(σ.p_c) decomposition
+                    // S: dot product
+                    ref[0] += c0123 * (D(X, X) + D(Y, Y) + D(Z, Z));
+                    // X: cross_x
+                    ref[1] += c0123 * (D(Y, Z) - D(Z, Y));
+                    // Y: cross_y
+                    ref[2] += c0123 * (D(Z, X) - D(X, Z));
+                    // Z: cross_z
+                    ref[3] += c0123 * (D(X, Y) - D(Y, X));
+                  }
+                }
+              }
+            }
+
+            const double ABSOLUTE_DEVIATION_THRESHOLD = 5.0E-14;
+            const double RELATIVE_DEVIATION_THRESHOLD = 1.0E-9;
+
+            for (auto comp = 0; comp < 4; ++comp) {
+              auto abs_err = abs(ref[comp] - results[comp][ijkl]);
+              auto rel_err = abs(abs_err / (ref[comp] + 1e-50));
+
+              bool not_ok = rel_err > RELATIVE_DEVIATION_THRESHOLD &&
+                            abs_err > ABSOLUTE_DEVIATION_THRESHOLD;
+
+              if (not_ok) {
+                std::cout << "(SL|SL) (" << s0 << " " << s1 << " | " << s2
+                          << " " << s3 << ") elem=" << ijkl << " comp=" << comp
+                          << " ref=" << ref[comp]
+                          << " sp=" << results[comp][ijkl]
+                          << " abs_err=" << abs_err << std::endl;
+              }
+              REQUIRE(!not_ok);
+            }
+
+            ++ijkl;
+            END_FOR_CART
+            END_FOR_CART
+            END_FOR_CART
+            END_FOR_CART
+          }
+        }
+      }
+    }
+  }
+
   SECTION("3-center RKB (K|LS)") {
     std::vector<Shell> dfbs{
         Shell{{3.0}, {{0, false, {1.0}}}, {{0.0, 0.0, 0.0}}},

include/libint2/sigma_p.h

@@ -51,9 +51,12 @@ inline SigmaPShell sigma_p(const Shell& s) { return {s}; }
 ///   BraKet::xs_xx — 3-center (K | a b)
 ///   BraKet::xx_xs — 3-center (a b | K)
 ///
-/// Component ordering: index = bra_q * n_ket_comp + ket_q
-///   For 2 (σ.p) on a particle: quaternion {S=0, X=1, Y=2, Z=3}
-///   For 1 (σ.p): Cartesian {x=0, y=1, z=2}
+/// Component structure based on total (σ.p) count:
+///   0 (σ.p) → 1 component (standard integral)
+///   1 (σ.p) → 3 Cartesian components {x=0, y=1, z=2} (3-center DF only)
+///   2 (σ.p) → 4 quaternion components {S=0, X=1, Y=2, Z=3}
+///              (same or different particles, e.g. (LL|SS) or (SL|SL))
+///   4 (σ.p) → 16 components, index = 4*bra_q + ket_q (SS|SS)
 class SigmaPEngine {
  public:
   using real_t = double;
@@ -164,17 +167,20 @@ class SigmaPEngine {
     return cart_idx(am + delta, nx, ny);
   }
 
-  /// Number of output components for a particle with n (σ.p) centers
-  static int particle_ncomp(int n_sp) {
-    switch (n_sp) {
+  /// Number of output components based on TOTAL (σ.p) count.
+  /// 0 → 1, 1 → 3 Cartesian, 2 → 4 quaternion, 4 → 16 (two pairs)
+  static int total_ncomp(int total_sp) {
+    switch (total_sp) {
       case 0:
         return 1;
       case 1:
         return 3;
       case 2:
         return 4;
+      case 4:
+        return 16;
       default:
-        assert(false);
+        assert(false && "unsupported (σ.p) count");
         return 0;
     }
   }
@@ -210,18 +216,16 @@ class SigmaPEngine {
     const int nc = ncart(lc), nd = ncart(ld);
     const int n_out = na * nb * nc * nd;
 
-    const int bra_sp = (sp_a ? 1 : 0) + (sp_b ? 1 : 0);
-    const int ket_sp = (sp_c ? 1 : 0) + (sp_d ? 1 : 0);
-    const int n_bra = particle_ncomp(bra_sp);
-    const int n_ket = particle_ncomp(ket_sp);
-    const int n_comp = n_bra * n_ket;
+    const int total_sp =
+        (sp_a ? 1 : 0) + (sp_b ? 1 : 0) + (sp_c ? 1 : 0) + (sp_d ? 1 : 0);
+    const int n_comp = total_ncomp(total_sp);
 
     buf_.assign(n_comp * n_out, 0.0);
     targets_.resize(n_comp);
     for (int q = 0; q < n_comp; ++q) targets_[q] = buf_.data() + q * n_out;
 
     // No (σ.p) → standard ERI
-    if (n_comp == 1) {
+    if (total_sp == 0) {
       const auto& r = coulomb_engine_.compute(a, b, c, d);
       if (r[0]) std::memcpy(buf_.data(), r[0], n_out * sizeof(real_t));
       return targets_;
@@ -300,73 +304,65 @@ class SigmaPEngine {
     }
 
     // ---- Contraction ----
-    for (int bq = 0; bq < n_bra; ++bq) {
-      for (int kq = 0; kq < n_ket; ++kq) {
-        int comp = bq * n_ket + kq;
-        real_t* out = buf_.data() + comp * n_out;
-
-        // T4 sign (only for 2+2 (σ.p) case)
-        int t4 = (bra_sp == 2 && ket_sp == 2 && bq > 0 && kq > 0) ? -1 : 1;
-
-        // Get direction pairs for bra and ket
-        DirPair bp[3], kp[3];
-        int nbp, nkp;
-        get_particle_pairs(sp_a, sp_b, bra_sp, bq, bp, nbp);
-        get_particle_pairs(sp_c, sp_d, ket_sp, kq, kp, nkp);
-
-        for (int ib = 0; ib < nbp; ++ib) {
-          // Directions: r[i] for center i. -1 means "not (σ.p), no direction"
-          int r[4];
-          r[0] = sp_a ? bp[ib].r1 : -1;
-          r[1] = sp_b ? bp[ib].r2 : -1;
-          const int bsign = bp[ib].sign;
-
-          for (int ik = 0; ik < nkp; ++ik) {
-            r[2] = sp_c ? kp[ik].r1 : -1;
-            r[3] = sp_d ? kp[ik].r2 : -1;
-            const int overall_sign = bsign * kp[ik].sign * t4;
-
-            // Loop over AM-shift combinations for each (σ.p) center
-            // delta[i]: +1 (up) or -1 (down). 0 if no (σ.p).
-            // di: 0=up, 1=down in the eri index
-            for (int d0 = 0; d0 <= (sp[0] ? 1 : 0); ++d0) {
-              if (d0 == 1 && !dn[0]) continue;
-              const int delta0 = sp[0] ? (d0 == 0 ? +1 : -1) : 0;
-
-              for (int d1 = 0; d1 <= (sp[1] ? 1 : 0); ++d1) {
-                if (d1 == 1 && !dn[1]) continue;
-                const int delta1 = sp[1] ? (d1 == 0 ? +1 : -1) : 0;
-
-                for (int d2 = 0; d2 <= (sp[2] ? 1 : 0); ++d2) {
-                  if (d2 == 1 && !dn[2]) continue;
-                  const int delta2 = sp[2] ? (d2 == 0 ? +1 : -1) : 0;
-
-                  for (int d3 = 0; d3 <= (sp[3] ? 1 : 0); ++d3) {
-                    if (d3 == 1 && !dn[3]) continue;
-                    const int delta3 = sp[3] ? (d3 == 0 ? +1 : -1) : 0;
-
-                    // Down-shifts contribute -1 sign each
-                    const int shift_sign =
-                        (delta0 < 0 ? -1 : 1) * (delta1 < 0 ? -1 : 1) *
-                        (delta2 < 0 ? -1 : 1) * (delta3 < 0 ? -1 : 1);
-
-                    const int eidx = d0 * 8 + d1 * 4 + d2 * 2 + d3;
-                    const real_t* e = eri[eidx];
-                    if (!e) continue;
-
-                    const int sns = eri_nb[eidx] * eri_nc[eidx] * eri_nd[eidx];
-                    const int snc = eri_nc[eidx] * eri_nd[eidx];
-                    const int snd = eri_nd[eidx];
-
-                    const int deltas[4] = {delta0, delta1, delta2, delta3};
-                    const int sign = overall_sign * shift_sign;
-
-                    // Inner loop over output CGFs
-                    accumulate(out, e, sign, shells, sp, am, r, deltas, na, nb,
-                               nc, nd, sns, snc, snd);
-                  }
-                }
-              }
+    // Identify which centers have (σ.p) — these form quaternion pairs.
+    // sp_centers[i] = index of i-th (σ.p) center in {0=a,1=b,2=c,3=d}
+    int sp_centers[4], n_sp_centers = 0;
+    for (int i = 0; i < 4; ++i)
+      if (sp[i]) sp_centers[n_sp_centers++] = i;
+
+    if (total_sp == 1) {
+      // Single unpaired (σ.p): 3 Cartesian components (DF only)
+      const int ci = sp_centers[0];  // the (σ.p) center
+      for (int r_dir = 0; r_dir < 3; ++r_dir) {
+        real_t* out = buf_.data() + r_dir * n_out;
+        int r[4] = {-1, -1, -1, -1};
+        r[ci] = r_dir;
+        contract_shifts(out, eri, eri_nb, eri_nc, eri_nd, sp, dn, am, r, 1, na,
+                        nb, nc, nd);
+      }
+    } else if (total_sp == 2) {
+      // One quaternion pair: 4 components {S, X, Y, Z}
+      const int c1 = sp_centers[0], c2 = sp_centers[1];
+      for (int q = 0; q < 4; ++q) {
+        real_t* out = buf_.data() + q * n_out;
+        DirPair pairs[3];
+        int npairs;
+        get_qpairs(q, pairs, npairs);
+        for (int ip = 0; ip < npairs; ++ip) {
+          int r[4] = {-1, -1, -1, -1};
+          r[c1] = pairs[ip].r1;
+          r[c2] = pairs[ip].r2;
+          contract_shifts(out, eri, eri_nb, eri_nc, eri_nd, sp, dn, am, r,
+                          pairs[ip].sign, na, nb, nc, nd);
+        }
+      }
+    } else if (total_sp == 4) {
+      // Two quaternion pairs: 16 components, index = 4*bra_q + ket_q
+      // Bra pair: first two (σ.p) centers, Ket pair: last two
+      // With standard ordering (a,b,c,d), bra=(0,1), ket=(2,3)
+      const int b1 = sp_centers[0], b2 = sp_centers[1];
+      const int k1 = sp_centers[2], k2 = sp_centers[3];
+      for (int bq = 0; bq < 4; ++bq) {
+        DirPair bp[3];
+        int nbp;
+        get_qpairs(bq, bp, nbp);
+        for (int kq = 0; kq < 4; ++kq) {
+          int comp = bq * 4 + kq;
+          real_t* out = buf_.data() + comp * n_out;
+          // T4 sign: -1 when both quaternion indices are cross-product
+          int t4 = (bq > 0 && kq > 0) ? -1 : 1;
+          DirPair kp[3];
+          int nkp;
+          get_qpairs(kq, kp, nkp);
+          for (int ib = 0; ib < nbp; ++ib) {
+            for (int ik = 0; ik < nkp; ++ik) {
+              int r[4] = {-1, -1, -1, -1};
+              r[b1] = bp[ib].r1;
+              r[b2] = bp[ib].r2;
+              r[k1] = kp[ik].r1;
+              r[k2] = kp[ik].r2;
+              contract_shifts(out, eri, eri_nb, eri_nc, eri_nd, sp, dn, am, r,
+                              bp[ib].sign * kp[ik].sign * t4, na, nb, nc, nd);
             }
           }
         }
@@ -376,30 +372,53 @@ class SigmaPEngine {
     return targets_;
   }
 
-  /// Get direction pairs for a particle based on (σ.p) config
-  static void get_particle_pairs(bool sp_f1, bool sp_f2, int n_sp, int q_idx,
-                                 DirPair* pairs, int& npairs) {
-    if (n_sp == 0) {
-      pairs[0] = {0, 0, 1};
-      npairs = 1;
-    } else if (n_sp == 1) {
-      // Single (σ.p): q_idx is the Cartesian direction (x=0, y=1, z=2)
-      if (sp_f1) {
-        pairs[0] = {q_idx, 0, 1};
-      } else {
-        pairs[0] = {0, q_idx, 1};
+  /// Contract AM-shift terms for a fixed set of directions r[0..3].
+  /// Loops over all up/down shift combinations and accumulates into out.
+  void contract_shifts(real_t* out, const std::array<const real_t*, 16>& eri,
+                       const std::array<int, 16>& eri_nb,
+                       const std::array<int, 16>& eri_nc,
+                       const std::array<int, 16>& eri_nd, const bool sp[4],
+                       const Shell* dn[4], const int am[4], const int r[4],
+                       int dir_sign, int na, int nb, int nc, int nd) {
+    for (int d0 = 0; d0 <= (sp[0] ? 1 : 0); ++d0) {
+      if (d0 == 1 && !dn[0]) continue;
+      const int delta0 = sp[0] ? (d0 == 0 ? +1 : -1) : 0;
+      for (int d1 = 0; d1 <= (sp[1] ? 1 : 0); ++d1) {
+        if (d1 == 1 && !dn[1]) continue;
+        const int delta1 = sp[1] ? (d1 == 0 ? +1 : -1) : 0;
+        for (int d2 = 0; d2 <= (sp[2] ? 1 : 0); ++d2) {
+          if (d2 == 1 && !dn[2]) continue;
+          const int delta2 = sp[2] ? (d2 == 0 ? +1 : -1) : 0;
+          for (int d3 = 0; d3 <= (sp[3] ? 1 : 0); ++d3) {
+            if (d3 == 1 && !dn[3]) continue;
+            const int delta3 = sp[3] ? (d3 == 0 ? +1 : -1) : 0;
+
+            const int shift_sign =
+                (delta0 < 0 ? -1 : 1) * (delta1 < 0 ? -1 : 1) *
+                (delta2 < 0 ? -1 : 1) * (delta3 < 0 ? -1 : 1);
+
+            const int eidx = d0 * 8 + d1 * 4 + d2 * 2 + d3;
+            const real_t* e = eri[eidx];
+            if (!e) continue;
+
+            const int sns = eri_nb[eidx] * eri_nc[eidx] * eri_nd[eidx];
+            const int snc = eri_nc[eidx] * eri_nd[eidx];
+            const int snd = eri_nd[eidx];
+            const int deltas[4] = {delta0, delta1, delta2, delta3};
+
+            accumulate(out, e, dir_sign * shift_sign, sp, am, r, deltas, na, nb,
+                       nc, nd, sns, snc, snd);
+          }
+        }
       }
-      npairs = 1;
-    } else {
-      get_qpairs(q_idx, pairs, npairs);
     }
   }
 
   /// Accumulate one AM-shift term into the output buffer
   void accumulate(real_t* out, const real_t* eri_buf, int sign,
-                  const Shell* shells[4], const bool sp[4], const int am[4],
-                  const int r[4], const int delta[4], int na, int nb, int nc,
-                  int nd, int sns, int snc, int snd) {
+                  const bool sp[4], const int am[4], const int r[4],
+                  const int delta[4], int na, int nb, int nc, int nd, int sns,
+                  int snc, int snd) {
     int oa = 0;
     for (int ax = am[0]; ax >= 0; --ax) {
       for (int ay = am[0] - ax; ay >= 0; --ay) {