Follow-ups to the reviews of #399: lmax checks, hartree-fock++ hessian coverage, compiler cleanups

.github/workflows/cmake.yml

@@ -45,15 +45,18 @@ jobs:
             build_shared_libs: ON
             os: macos-15
             do_fortran: ON
-          # exercise a center-dependent AM table for the 3-center RKB ERIs and
-          # 2nd geometric derivatives of the 4-center ERIs in one lane only
+          # exercise a center-dependent AM table for the 3-center RKB ERIs,
+          # per-derivative AM lists, and 2nd geometric derivatives of the
+          # 4-center ERIs in one lane only; the 2nd-derivative ERIs must cover
+          # the d shells of the basis of hartree-fock++, else its hessian is
+          # skipped
           - build_type: Release
             os: macos-15
             subproject: OFF
             rkb_eri3_max_am: 3
             eri_deriv_order: 2
             max_am: '2;2;1'
-            eri_max_am: '2;2;1'
+            eri_max_am: '2;2;2'
 
     name: "Repo • ${{ matrix.os }}: ${{ matrix.cxx }} ${{ matrix.build_type }} shared=${{ matrix.build_shared_libs }} fetch=${{ matrix.subproject }}"
     runs-on: ${{ matrix.os }}

CHANGES

@@ -4,6 +4,17 @@
 Following is a brief summary of changes made in each release of Libint.
 
 - unreleased:
+  - PR #450: Engine::compute() throws Engine::lmax_exceeded, in all builds, if a shell exceeds the
+             max_l of the Engine (previously a Debug-only assert); Engine::lmax_exceeded::lmax_limit()
+             is the max angular momentum also when thrown by the constructor/set() (previously the
+             max angular momentum + 1)
+  - PR #450: Engine::set(BraKet) and Engine::set(Operator) resize the primitive data for the rank
+             of the new braket (previously e.g. switching an Engine from a 3-center to a 4-center
+             braket overran a buffer); Engine::max_nprim() returns the max number of primitives the
+             Engine was sized for (previously derived from the last computed shell pair)
+  - PR #450: the generated Operator::coulomb_op kernels compute each component in place, without
+             an element-wise copy of the derivative integrals (the compiler makes the target of a
+             recurrence relation that is an alias of its child refer to the child)
   - PR #399: relativistic restricted kinetic balance (RKB) integrals: 2-body Operator::coulomb_opop,
     Operator::coulomb_op, Operator::opop_coulomb_opop, Operator::op_coulomb_op (enabled by
     LIBINT2_ENABLE_RKB_ERI and LIBINT2_ENABLE_RKB_ERI3, off by default; no geometric derivatives yet)

cmake/modules/int_am.cmake

@@ -357,6 +357,15 @@ endmacro()
 
 process_integrals_class(ONEBODY)
 process_integrals_class(ERI)
+# max AM of the 4-center ERIs for each derivative order, comma-separated
+# (for the features file read by the Hartree-Fock tests)
+set(Libint2_MAX_AM_ERI_DERIVS "")
+if (LIBINT_INCLUDE_ERI GREATER -1)
+    foreach(_d RANGE ${LIBINT_INCLUDE_ERI})
+        list(APPEND Libint2_MAX_AM_ERI_DERIVS ${_candidate_ERI_d${_d}})
+    endforeach()
+endif()
+list(JOIN Libint2_MAX_AM_ERI_DERIVS "," Libint2_MAX_AM_ERI_DERIVS)
 process_integrals_class(RKB_ERI)
 process_integrals_class(RKB_ERI3)
 process_integrals_class(ERI3)

cmake/modules/int_computed.cmake.in

@@ -87,3 +87,4 @@ list(PREPEND Libint2_CONFIG_COMPONENTS "${_o1st}${_o2nd}")
 
 # for features and libint2-config.cmake files
 set(Libint2_MAX_AM_ERI "@_candidate_ERI_d0@")
+set(Libint2_MAX_AM_ERI_DERIVS "@Libint2_MAX_AM_ERI_DERIVS@")

export/fortran/test-eri.cc

@@ -5,6 +5,8 @@
 #include <libint2/util/generated/libint2_types.h>
 #include <libint2/util/memory.h>
 
+#include <algorithm>
+#include <cmath>
 #include <iomanip>
 
 #include "../tests/eri/prep_libint2.h"
@@ -498,13 +500,20 @@ bool test_eri_c_f(int &contrdepth, int &am1, double *c1, double *alpha1,
     const unsigned int n2 = (am2 + 1) * (am2 + 2) / 2;
     const unsigned int n3 = (am3 + 1) * (am3 + 2) / 2;
     const unsigned int n4 = (am4 + 1) * (am4 + 2) / 2;
+    // the C and Fortran codes compute the primitive data with differently
+    // ordered arithmetic; the resulting roundoff is relative to the magnitude
+    // of the shell set, not of each element (elements of derivative shell sets
+    // can be small differences of large terms)
+    double max_abs = 0.0;
+    for (unsigned int i = 0; i != n1 * n2 * n3 * n4; ++i)
+      max_abs = std::max(max_abs, std::abs(eri_shell_set_c[i]));
     unsigned int nel = 0;
     for (int a = 0; a < n1; a++) {
       for (int b = 0; b < n2; b++) {
         for (int c = 0; c < n3; c++) {
           for (int d = 0; d < n4; d++, nel++) {
             const double abs_error = abs(*eri_shell_set_c - *eri_shell_set_f);
-            const double rel_error = abs_error / abs(*eri_shell_set_c);
+            const double rel_error = abs_error / max_abs;
             if (rel_error > threshold) {
               std::cout << std::setprecision(17) << "Elem " << nel
                         << " di= " << deriv_order

export/tests/hartree-fock/hartree-fock++-validate.py

@@ -44,6 +44,21 @@ def validate(label, data, refdata, tolerance, textline):
 if LIBINT_ERI_MAX_AM<2:
     sys.exit(0)
 
+# the max AM of the basis (aug-cc-pVDZ for H2O)
+OBS_MAX_L = 2
+
+def skip_ok(label, deriv_order, textline):
+    """a quantity that needs the ERI derivatives of order deriv_order may only
+    be skipped if their kernels do not support the AM of the basis"""
+    expected = LIBINT_ERI_MAX_AM_DERIVS[deriv_order] < OBS_MAX_L
+    if expected:
+        print(textline, end="")
+    else:
+        print(label, "check: failed\nskipped although the library supports l <=",
+              LIBINT_ERI_MAX_AM_DERIVS[deriv_order], "for derivative order",
+              deriv_order, "\nactual:", textline)
+    return expected
+
 eref = [-76.003354058439]
 etol = 5e-12
 
@@ -175,17 +190,13 @@ def validate(label, data, refdata, tolerance, textline):
     elif match15:
         sQok = validate("spherical electric quadrupole moment", match15.groups(), sQref, sQtol, line)
     elif match16:
-        F2ok = True
-        print(line,end="")
+        F2ok = skip_ok("2-body force", 1, line)
     elif match17:
-        Fok = True
-        print(line,end="")
+        Fok = skip_ok("HF force", 1, line)
     elif match18:
-        H2ok = True
-        print(line,end="")
+        H2ok = skip_ok("2-body hessian", 2, line)
     elif match19:
-        Hok = True
-        print(line,end="")
+        Hok = skip_ok("HF hessian", 2, line)
     else:
         print(line,end="")
 

export/tests/hartree-fock/hartree-fock++.cc

@@ -2589,8 +2589,14 @@ void api_basic_compile_test(const BasisSet& obs,
   }
 
 #if LIBINT2_DERIV_ERI_ORDER
-  // skipped if the derivative kernels do not support the AM of the basis
-  if (obs.max_l() <= LIBINT2_MAX_AM_eri1) {  // test deriv 2-index ints
+  // skipped if the derivative kernels do not support the AM of the basis; this
+  // needs the 4-center (eri1) and the 2-center (2eri1) kernels
+#if defined(LIBINT2_MAX_AM_2eri1)
+  const int lmax_d1 = std::min(LIBINT2_MAX_AM_eri1, LIBINT2_MAX_AM_2eri1);
+#else
+  const int lmax_d1 = LIBINT2_MAX_AM_eri1;
+#endif
+  if (obs.max_l() <= lmax_d1) {  // test deriv 2-index ints
     Engine eri4_engine(Operator::coulomb, obs.max_nprim(), obs.max_l(), 1);
     Engine eri2_engine = eri4_engine;
     eri2_engine.set(BraKet::xs_xs);
@@ -2625,7 +2631,12 @@ void api_basic_compile_test(const BasisSet& obs,
 #endif
 
 #if LIBINT2_DERIV_ERI_ORDER > 1
-  if (obs.max_l() <= LIBINT2_MAX_AM_eri2) {  // test 2nd deriv 2-index ints
+#if defined(LIBINT2_MAX_AM_2eri2)
+  const int lmax_d2 = std::min(LIBINT2_MAX_AM_eri2, LIBINT2_MAX_AM_2eri2);
+#else
+  const int lmax_d2 = LIBINT2_MAX_AM_eri2;
+#endif
+  if (obs.max_l() <= lmax_d2) {  // test 2nd deriv 2-index ints
     Engine eri4_engine(Operator::coulomb, obs.max_nprim(), obs.max_l(), 2);
     Engine eri2_engine = eri4_engine;
     eri2_engine.set(BraKet::xs_xs);

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

@@ -21,6 +21,7 @@
 #include <libint2/deriv_iter.h>
 
 #include <algorithm>
+#include <functional>
 #include <limits>
 #include <sstream>
 
@@ -741,6 +742,82 @@ TEST_CASE("rejected set() leaves no kernels", "[engine][2-body]") {
 #endif
 }
 
+// compute() must throw Engine::lmax_exceeded, in all builds, for a shell
+// beyond the max_l of the Engine or beyond the limit of the kernels for its
+// center
+TEST_CASE("compute() checks the angular momentum limits", "[engine][2-body]") {
+  const auto eps = std::numeric_limits<double>::epsilon();
+  const Shell s{{1.0}, {{0, false, {1.0}}}, {{0.0, 0.0, 0.0}}};
+  // @return the lmax_exceeded thrown by @p f
+  auto thrown_by =
+      [](const std::function<void()> &f) -> std::pair<size_t, size_t> {
+    try {
+      f();
+    } catch (Engine::lmax_exceeded &e) {
+      return {e.lmax_limit(), e.lmax_requested()};
+    }
+    FAIL("lmax_exceeded was not thrown");
+    return {0, 0};
+  };
+
+#if defined(LIBINT2_SUPPORT_ERI) && LIBINT2_MAX_AM_eri > 0
+  SECTION("max_l of the Engine") {
+    const Shell p{{1.0}, {{1, false, {1.0}}}, {{0.0, 0.0, 1.0}}};
+    Engine engine(Operator::coulomb, 1, 0);
+    CHECK(thrown_by([&]() { engine.compute(s, s, s, p); }) ==
+          std::make_pair(size_t(0), size_t(1)));
+    Engine engine1(Operator::overlap, 1, 0);
+    CHECK(thrown_by([&]() { engine1.compute(p, s); }) ==
+          std::make_pair(size_t(0), size_t(1)));
+  }
+#endif
+
+  // the kernels of a center-dependent task support higher AM on the fitting
+  // center than on the AO centers, the Engine's max_l does not tell them apart
+#if defined(LIBINT2_SUPPORT_ERI3) && LIBINT2_CENTER_DEPENDENT_MAX_AM_3eri && \
+    LIBINT2_MAX_AM_3eri > LIBINT2_MAX_AM_default
+  SECTION("AM limit of the AO centers of 3-center kernels") {
+    const int l_ao = LIBINT2_MAX_AM_default + 1;
+    const Shell ao{{1.0}, {{l_ao, true, {1.0}}}, {{0.0, 0.0, 1.0}}};
+    Engine engine(Operator::coulomb, 1, LIBINT2_MAX_AM_3eri, 0, eps,
+                  libint2::default_params(Operator::coulomb), BraKet::xs_xx);
+    CHECK(thrown_by([&]() { engine.compute(s, ao, s); }) ==
+          std::make_pair(size_t(LIBINT2_MAX_AM_default), size_t(l_ao)));
+  }
+#endif
+
+  // the limit reported by initialize() is the max AM, not the max AM + 1
+#if defined(LIBINT2_SUPPORT_ERI)
+  SECTION("AM limit reported by initialize()") {
+    CHECK(thrown_by([&]() {
+            Engine engine(Operator::coulomb, 1, LIBINT2_MAX_AM_eri + 1);
+          }) == std::make_pair(size_t(LIBINT2_MAX_AM_eri),
+                               size_t(LIBINT2_MAX_AM_eri + 1)));
+  }
+#endif
+}
+
+// set(BraKet) must resize the core evaluator for the new braket: a 3-center
+// Engine needs Boys functions up to 3·max_l, a 4-center one up to 4·max_l
+TEST_CASE("set(BraKet) rebuilds the core evaluator", "[engine][2-body]") {
+#if defined(LIBINT2_SUPPORT_ERI) && defined(LIBINT2_SUPPORT_ERI3)
+  const int l = std::min(LIBINT2_MAX_AM_eri, LIBINT2_MAX_AM_3eri);
+  const Shell a{{1.0, 0.3}, {{l, false, {0.4, 0.7}}}, {{0.0, 0.0, 0.0}}};
+  const Shell b{{1.1, 0.2}, {{l, false, {0.5, 0.6}}}, {{0.0, 0.5, 0.1}}};
+  const auto eps = std::numeric_limits<double>::epsilon();
+  Engine engine(Operator::coulomb, 2, l, 0, eps,
+                libint2::default_params(Operator::coulomb), BraKet::xs_xx);
+  engine.set(BraKet::xx_xx);
+  Engine ref(Operator::coulomb, 2, l);
+  const auto &result = engine.compute(a, b, b, a);
+  const auto &result_ref = ref.compute(a, b, b, a);
+  REQUIRE(result[0] != nullptr);
+  const auto n = a.size() * b.size() * b.size() * a.size();
+  for (size_t i = 0; i != n; ++i)
+    REQUIRE(result[0][i] == Approx(result_ref[0][i]).margin(1e-14));
+#endif
+}
+
 TEST_CASE("RKB Coulomb integrals", "[engine][2-body]") {
   const std::vector<Shell> all_shells{
       // pseudorandom s

include/libint2/engine.h

@@ -54,6 +54,7 @@
 #include <limits>
 #include <map>
 #include <memory>
+#include <string>
 #include <tuple>
 #include <utility>
 #include <vector>
@@ -308,8 +309,13 @@ enum class Operator {
  σ_i σ_j = δ_ij + i ε_ijk σ_k into a scalar (component 0) and a vector
  (components 1-3) part; the user assembles the 2×2 spin block of that electron
  as q₀·1 + i (q₁ σ_x + q₂ σ_y + q₃ σ_z) and must not apply the Pauli matrices
- again. "Not folded" means the result is the spatial tensor and the user
- contracts it with the Pauli matrices.
+ again. This applies as stated to the operators with one folded electron
+ (`coulomb_opop` and the 1-body operators); for `opop_coulomb_opop`, which
+ folds both electrons, use its prefactor column instead: the product of the two
+ i's of the vector ⊗ vector block (i² = −1) is already in the kernel, so
+ assembling each electron's block as q₀ + i σ·q would apply it twice.
+ "Not folded" means the result is the spatial tensor and the user contracts it
+ with the Pauli matrices.
 
  | Engine (C++ API)    | Compiler          | task label          | #comp | σ·p acts on | component index                                                                                                                                              | prefactor                                                                                                                                                                                                   |
  |---------------------|-------------------|---------------------|-------|-------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
@@ -903,6 +909,7 @@ class Engine {
         primdata_(std::move(other.primdata_)),
         spbra_(std::move(other.spbra_)),
         spket_(std::move(other.spket_)),
+        max_nprim_(other.max_nprim_),
         stack_size_(other.stack_size_),
         lmax_(other.lmax_),
         hard_lmax_(other.hard_lmax_),
@@ -935,6 +942,7 @@ class Engine {
         primdata_(other.primdata_.size()),
         spbra_(other.spbra_),
         spket_(other.spket_),
+        max_nprim_(other.max_nprim_),
         stack_size_(other.stack_size_),
         lmax_(other.lmax_),
         deriv_order_(other.deriv_order_),
@@ -958,6 +966,7 @@ class Engine {
     primdata_ = std::move(other.primdata_);
     spbra_ = std::move(other.spbra_);
     spket_ = std::move(other.spket_);
+    max_nprim_ = other.max_nprim_;
     stack_size_ = other.stack_size_;
     lmax_ = other.lmax_;
     hard_lmax_ = other.hard_lmax_;
@@ -991,6 +1000,7 @@ class Engine {
     primdata_.resize(other.primdata_.size());
     spbra_ = other.spbra_;
     spket_ = other.spket_;
+    max_nprim_ = other.max_nprim_;
     stack_size_ = other.stack_size_;
     lmax_ = other.lmax_;
     deriv_order_ = other.deriv_order_;
@@ -1069,17 +1079,14 @@ class Engine {
   }
 
   /// @return the maximum number of primitives that this engine can handle
-  std::size_t max_nprim() const {
-    assert(spbra_.primpairs.size() == spket_.primpairs.size());
-    return static_cast<std::size_t>(std::sqrt(spbra_.primpairs.size()));
-  }
+  std::size_t max_nprim() const { return max_nprim_; }
 
   /// reset the maximum number of primitives
   /// @param[in] n the maximum number of primitives
   /// @note left unchanged if the value returned by Engine::max_nprim is greater
   /// than @c n
   Engine& set_max_nprim(std::size_t n) {
-    if (n * n > spbra_.primpairs.size()) {
+    if (n > max_nprim_) {
       spbra_.resize(n);
       spket_.resize(n);
       initialize(n);
@@ -1273,10 +1280,18 @@ class Engine {
   /// Exception class to be used when the angular momentum limit is exceeded.
   class lmax_exceeded : virtual public std::logic_error {
    public:
+    /// @param task_name the label of the task (or the empty string)
+    /// @param lmax_limit the max angular momentum supported
+    /// @param lmax_requested the angular momentum that was requested
     lmax_exceeded(const char* task_name, size_t lmax_limit,
                   size_t lmax_requested)
         : std::logic_error(
-              "Engine::lmax_exceeded -- angular momentum limit exceeded"),
+              std::string(
+                  "Engine::lmax_exceeded -- angular momentum limit exceeded") +
+              (task_name[0] != '\0' ? std::string(" for task ") + task_name
+                                    : std::string()) +
+              ": the max angular momentum is " + std::to_string(lmax_limit) +
+              ", requested " + std::to_string(lmax_requested)),
           lmax_limit_(lmax_limit),
           lmax_requested_(lmax_requested) {
       strncpy(task_name_, task_name, 64);
@@ -1324,6 +1339,8 @@ class Engine {
   BraKet braket_;
   std::vector<Libint_t> primdata_;
   ShellPair spbra_, spket_;
+  // primdata_ is sized for max_nprim_ primitives per shell; set by initialize()
+  std::size_t max_nprim_ = 0;
   size_t stack_size_;  // amount allocated by libint2_init_xxx in
                        // primdata_[0].stack
   int lmax_;