diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml
index e10e10a3c..5cf60d69a 100644
--- a/.github/workflows/cmake.yml
+++ b/.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 }}
diff --git a/CHANGES b/CHANGES
index e76e390c8..c48596a6f 100644
--- a/CHANGES
+++ b/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)
diff --git a/cmake/modules/int_am.cmake b/cmake/modules/int_am.cmake
index cabcc18bd..4f0ee80e9 100644
--- a/cmake/modules/int_am.cmake
+++ b/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)
diff --git a/cmake/modules/int_computed.cmake.in b/cmake/modules/int_computed.cmake.in
index 51d79d527..4a66b89e2 100644
--- a/cmake/modules/int_computed.cmake.in
+++ b/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@")
diff --git a/export/fortran/test-eri.cc b/export/fortran/test-eri.cc
index 50cc86306..7cf32dad8 100644
--- a/export/fortran/test-eri.cc
+++ b/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
diff --git a/export/tests/hartree-fock/hartree-fock++-validate.py b/export/tests/hartree-fock/hartree-fock++-validate.py
index 95be09631..8e6f2a929 100644
--- a/export/tests/hartree-fock/hartree-fock++-validate.py
+++ b/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="")
 
diff --git a/export/tests/hartree-fock/hartree-fock++.cc b/export/tests/hartree-fock/hartree-fock++.cc
index 2e45d3b41..ac45d1d9c 100644
--- a/export/tests/hartree-fock/hartree-fock++.cc
+++ b/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);
diff --git a/export/tests/unit/test-2body.cc b/export/tests/unit/test-2body.cc
index d346135e1..621648339 100644
--- a/export/tests/unit/test-2body.cc
+++ b/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
diff --git a/include/libint2/engine.h b/include/libint2/engine.h
index badd60654..1d3d79cc0 100644
--- a/include/libint2/engine.h
+++ b/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_;
diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h
index e1a90f254..1b254eaaa 100644
--- a/include/libint2/engine.impl.h
+++ b/include/libint2/engine.impl.h
@@ -21,6 +21,7 @@
 #ifndef _libint2_src_lib_libint_engineimpl_h_
 #define _libint2_src_lib_libint_engineimpl_h_
 
+#include <initializer_list>
 #include <iterator>
 
 #include "./deriv_map.h"
@@ -217,8 +218,9 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute1(
 
   const auto l1 = s1.contr[0].l;
   const auto l2 = s2.contr[0].l;
-  assert(l1 <= lmax_ && "the angular momentum limit is exceeded");
-  assert(l2 <= lmax_ && "the angular momentum limit is exceeded");
+  // the scratch and the stack are sized for lmax_, check in all builds
+  if (std::max(l1, l2) > lmax_)
+    throw lmax_exceeded(task_name_, lmax_, std::max(l1, l2));
 
   // if want nuclear, make sure there is at least one nucleus .. otherwise the
   // user likely forgot to call set_params
@@ -274,7 +276,6 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute1(
 
   // adjust max angular momentum, if needed
   const auto lmax = std::max(l1, l2);
-  assert(lmax <= lmax_ && "the angular momentum limit is exceeded");
 
   // N.B. for l=0 no need to transform to solid harmonics
   // this is a workaround for the corner case of oper_ == Operator::*nuclear,
@@ -660,7 +661,7 @@ __libint2_engine_inline void Engine::_initialize() {
         std::max(hard_lmax_, hard_default_lmax_), hard_lmax_);                 \
     if (lmax_ >= lmax) {                                                       \
       buildfnptrs_ = nullptr; /* leave no table for kernels we cannot use */   \
-      throw Engine::lmax_exceeded(BOOST_PP_STRINGIZE(BOOST_PP_NBODYENGINE_MCR3_TASK(product)), lmax, lmax_);              \
+      throw Engine::lmax_exceeded(BOOST_PP_STRINGIZE(BOOST_PP_NBODYENGINE_MCR3_TASK(product)), lmax - 1, lmax_);          \
     }                                                                          \
     if (stack_size_ > 0)                                                       \
       LIBINT2_PREFIXED_NAME(libint2_cleanup_default)(&primdata_[0]);           \
@@ -720,7 +721,10 @@ __libint2_engine_inline void Engine::initialize(size_t max_nprim) {
   // initialize braket, if needed
   if (braket_ == BraKet::invalid) braket_ = default_braket(oper_);
 
-  if (max_nprim != 0) primdata_.resize(std::pow(max_nprim, braket_rank()));
+  // the size of primdata_ depends on the rank of the braket, hence resize it
+  // also when called without max_nprim, e.g. by set(BraKet) or set(Operator)
+  if (max_nprim != 0) max_nprim_ = max_nprim;
+  primdata_.resize(std::pow(max_nprim_, braket_rank()));
 
   // initialize targets
   {
@@ -1336,6 +1340,10 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
           (tspbra != nullptr && tspket != nullptr)) &&
          "Engine::compute2 -- expects zero or two ShellPair objects");
   assert(screening_method_ != ScreeningMethod::Invalid);
+  // the derivative order of the integrals over the operator, known at compile
+  // time, hence evaluated once rather than per primitive
+  constexpr auto total_deriv_order =
+      der + operator_traits<op>::intrinsic_deriv_order;
   if (buildfnptrs_ == nullptr) throw missing_kernel();
 
   //
@@ -1347,13 +1355,12 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
           tket2.ncontr() == 1) &&
          "generally-contracted shells are not yet supported");
 
-  // angular momentum limit obeyed? can only be fully checked in
-  // braket-dependent code, here only do a basic test that does not guarantee
-  // that the shell-set can be computed
-  assert(tbra1.contr[0].l <= lmax_ && "the angular momentum limit is exceeded");
-  assert(tbra2.contr[0].l <= lmax_ && "the angular momentum limit is exceeded");
-  assert(tket1.contr[0].l <= lmax_ && "the angular momentum limit is exceeded");
-  assert(tket2.contr[0].l <= lmax_ && "the angular momentum limit is exceeded");
+  // angular momentum limit obeyed? the scratch, the stack and the core
+  // evaluator are sized for lmax_, hence check in all builds; the limits of the
+  // kernels can only be fully checked in braket-dependent code
+  for (const auto* sh : {&tbra1, &tbra2, &tket1, &tket2})
+    if (sh->contr[0].l > lmax_)
+      throw lmax_exceeded(task_name_, lmax_, sh->contr[0].l);
 
 #if LIBINT2_SHELLQUARTET_SET == \
     LIBINT2_SHELLQUARTET_SET_STANDARD  // standard angular momentum ordering
@@ -1437,7 +1444,6 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
   // adjust max angular momentum, if needed
   auto lmax = std::max(std::max(bra1.contr[0].l, bra2.contr[0].l),
                        std::max(ket1.contr[0].l, ket2.contr[0].l));
-  assert(lmax <= lmax_ && "the angular momentum limit is exceeded");
   const auto lmax_bra = std::max(bra1.contr[0].l, bra2.contr[0].l);
   const auto lmax_ket = std::max(ket1.contr[0].l, ket2.contr[0].l);
 
@@ -1597,7 +1603,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_ + intrinsic_deriv_order();
+            const auto mmax = l + total_deriv_order;
 
             if (!skip_core_ints) {
 // Simple core-eval dispatch: just `core_eval_ptr->eval(gm_ptr, T, mmax)`.
@@ -1847,7 +1853,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
               const auto Wz =
                   (gammap_o_gammapgammaq * P[2] + gammaq_o_gammapgammaq * Q[2]);
 
-              if (deriv_order_ + intrinsic_deriv_order() > 0 || lmax_bra > 0) {
+              if (total_deriv_order > 0 || lmax_bra > 0) {
 #if LIBINT2_DEFINED(eri, WP_x)
                 primdata.WP_x[0] = Wx - P[0];
 #endif
@@ -1858,7 +1864,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
                 primdata.WP_z[0] = Wz - P[2];
 #endif
               }
-              if (deriv_order_ + intrinsic_deriv_order() > 0 || lmax_ket > 0) {
+              if (total_deriv_order > 0 || lmax_ket > 0) {
 #if LIBINT2_DEFINED(eri, WQ_x)
                 primdata.WQ_x[0] = Wx - Q[0];
 #endif
@@ -1942,7 +1948,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
 #endif
 
               // prefactors for derivative ERI relations
-              if (deriv_order_ + intrinsic_deriv_order() > 0) {
+              if (total_deriv_order > 0) {
 #if LIBINT2_DEFINED(eri, alpha1_rho_over_zeta2)
                 primdata.alpha1_rho_over_zeta2[0] =
                     alpha0 * (oogammap * gammaq_o_gammapgammaq);
@@ -2025,8 +2031,7 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
   }
 
   // compute directly (ss|ss)
-  const auto compute_directly =
-      lmax == 0 && deriv_order_ == 0 && intrinsic_deriv_order() == 0;
+  const auto compute_directly = lmax == 0 && total_deriv_order == 0;
 
   if (compute_directly) {
 #ifdef LIBINT2_ENGINE_TIMERS
diff --git a/src/bin/libint/build_libint.cc b/src/bin/libint/build_libint.cc
index 2726ea585..3f95d4c5e 100644
--- a/src/bin/libint/build_libint.cc
+++ b/src/bin/libint/build_libint.cc
@@ -31,9 +31,7 @@
   */
 
 #include <boost/preprocessor.hpp>
-#include <chrono>
 #include <fstream>
-#include <iomanip>
 #include <iostream>
 #include <limits>
 #include <vector>
@@ -95,7 +93,6 @@ struct ShellQuartetSetPredicate<ShellSetType_Standard> {
       case PairSwap::coupled:
         // the ket is sorted and the bra follows; when the totals tie the
         // engine swaps bra and ket such that max(la,lb) <= lc
-        assert(p1p2_swappable);
         return lc >= ld && (la + lb < lc + ld ||
                             (la + lb == lc + ld && std::max(la, lb) <= lc));
       case PairSwap::none:
@@ -113,7 +110,6 @@ struct ShellQuartetSetPredicate<ShellSetType_ORCA> {
         return la <= lb && lc <= ld &&
                (!p1p2_swappable || (la < lc || (la == lc && lb <= ld)));
       case PairSwap::coupled:
-        assert(p1p2_swappable);
         return lc <= ld && (la + lb > lc + ld ||
                             (la + lb == lc + ld && std::min(la, lb) >= lc));
       case PairSwap::none:
@@ -1312,6 +1308,10 @@ static void build_TwoPRep_2b_2k(
   // (see the descriptor's p1p2_swappable and pair_swap)
   constexpr bool p1p2_swappable = OperDescrType::p1p2_swappable;
   constexpr PairSwap pair_swap = OperDescrType::pair_swap;
+  // the shell-set predicates assume that coupled pair swaps come with a
+  // particle swap (the engine swaps the pairs together with bra and ket)
+  static_assert(pair_swap != PairSwap::coupled || p1p2_swappable,
+                "PairSwap::coupled requires a p1p2-swappable operator");
 
   // Note: la, lb, lc, ld generate code for chemist notation (ab|O|cd), where O
   // is a two-body operator.
@@ -1377,7 +1377,7 @@ static void build_TwoPRep_2b_2k(
           ////////////
           // loop over unique derivative index combinations
           ////////////
-          // NB translational invariance is now handled by CR_DerivGauss
+          // NB translational invariance is handled by DerivGaussV2
           CartesianDerivIterator<4> diter(deriv_level);
 
           bool last_deriv = false;
@@ -1609,7 +1609,7 @@ static void build_TwoPRep_1b_2k(
         ////////////
         // loop over unique derivative index combinations
         ////////////
-        // NB translational invariance is now handled by CR_DerivGauss
+        // NB translational invariance is handled by DerivGaussV2
         CartesianDerivIterator<3> diter(deriv_level);
         std::vector<std::shared_ptr<TwoBody_sh_11_11>> targets;
         bool last_deriv = false;
@@ -1832,7 +1832,7 @@ void build_TwoPRep_1b_1k(std::ostream& os,
       ////////////
       // loop over unique derivative index combinations
       ////////////
-      // NB translational invariance is now handled by CR_DerivGauss
+      // NB translational invariance is handled by DerivGaussV2
       CartesianDerivIterator<2> diter(deriv_level);
       std::vector<std::shared_ptr<TwoPRep_sh_11_11>> targets;
       bool last_deriv = false;
@@ -1909,7 +1909,6 @@ void build_TwoPRep_1b_1k(std::ostream& os,
 
       g_progress.current_task = label;
       g_progress.print();
-      std::cout.flush();
 
       std::string prefix(cparams->source_directory());
       std::deque<std::string> decl_filenames;
diff --git "a/src/bin/libint/comp_11_Coulomb\317\203p_11.h" "b/src/bin/libint/comp_11_Coulomb\317\203p_11.h"
index fd85b01f4..a0450eb36 100644
--- "a/src/bin/libint/comp_11_Coulomb\317\203p_11.h"
+++ "b/src/bin/libint/comp_11_Coulomb\317\203p_11.h"
@@ -25,11 +25,10 @@
 #include <entity.h>
 #include <gaussoper.h>
 #include <generic_rr.h>
+#include <rkb_fold_codegen.h>
 #include <task.h>
 #include <twoprep_11_11.h>
 
-#include <rkb_fold_codegen.h>
-
 #include <stdexcept>
 
 namespace libint2 {
@@ -64,6 +63,14 @@ class CR_11_Coulombσp_11 : public GenericRecurrenceRelation<
 
   static bool directional() { return false; }
 
+  /// each Cartesian component is a single derivative integral with unit
+  /// coefficient, hence the target is the child (see fold_terms())
+  bool is_alias() const override {
+    if (this->num_children() != 1) return false;
+    const auto terms = fold_terms(target_->oper()->descr().cartesian_index());
+    return terms.size() == 1 && terms[0].coef == 1;
+  }
+
  private:
   using ParentType::is_simple;
   using ParentType::target_;
diff --git "a/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h" "b/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h"
index 0404da395..212d2d164 100644
--- "a/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h"
+++ "b/src/bin/libint/comp_11_\317\203p\317\203pCoulomb\317\203p\317\203p_11.h"
@@ -25,11 +25,10 @@
 #include <entity.h>
 #include <gaussoper.h>
 #include <generic_rr.h>
+#include <rkb_fold_codegen.h>
 #include <task.h>
 #include <twoprep_11_11.h>
 
-#include <rkb_fold_codegen.h>
-
 #include <stdexcept>
 
 namespace libint2 {
@@ -48,8 +47,9 @@ namespace libint2 {
  * Molecules", Springer, 2007):
  *   \f$ (\sigma \cdot u)(\sigma \cdot v) = (u \cdot v)I
  *       + i\sigma \cdot (u \times v) \f$
- * applied independently to each electron's spin space gives a tensor product
- * of two quaternions with \f$ 4 \times 4 = 16 \f$ components:
+ * applied independently to each electron's spin space gives the tensor
+ * product of two 4-component (scalar + vector) objects, with
+ * \f$ 4 \times 4 = 16 \f$ components:
  *
  *   index = 4 * q₁ + q₂
  *
diff --git a/src/bin/libint/comp_deriv_gauss_v2.h b/src/bin/libint/comp_deriv_gauss_v2.h
index 172540677..018158826 100644
--- a/src/bin/libint/comp_deriv_gauss_v2.h
+++ b/src/bin/libint/comp_deriv_gauss_v2.h
@@ -436,8 +436,8 @@ std::string DerivGaussV2<IntType, part, where, trans_inv_part,
   if constexpr (using_trans_inv) {
     typedef typename TargetType::AuxIndexType mType;
     static std::shared_ptr<mType> aux0(new mType(0u));
-    // N.B. not "CR_DerivGauss...": CR_DerivGauss labels its own instances
-    // that way, and the RRStack deduplicates by label
+    // N.B. a prefix of its own, distinct from the "CR_DerivGauss " labels of
+    // CR_DerivGauss instances (the RRStack deduplicates by label)
     os << "DerivGaussV2TI"
        << "P" << part << to_string(where) << "D" << dir_
        << genintegralset_label(target_->bra(), target_->ket(), aux0,
diff --git a/src/bin/libint/dg.cc b/src/bin/libint/dg.cc
index 0955cbb05..3d9b901ac 100644
--- a/src/bin/libint/dg.cc
+++ b/src/bin/libint/dg.cc
@@ -846,6 +846,32 @@ void DirectedGraph::handle_trivial_nodes(
       }
     }
 
+    // Is the exit arc an RR whose target is identical to its child (e.g.
+    // CR_11_Coulombσp_11)? Then make the vertex refer to the child, so that the
+    // child is computed in its place and no code is generated for the RR.
+    // N.B. the child must have no other parent: a referring vertex is never
+    // allocated, hence never tags its child, and freeing it (once its parents
+    // are allocated) frees the child; if it is a target, the child takes over
+    // its address.
+    if (!(vptr)->refers_to_another()) {
+      std::shared_ptr<DGArcRR> arc_cast =
+          std::dynamic_pointer_cast<DGArcRR, DGArc>(arc);
+      if (arc_cast && arc_cast->rr()->is_alias()) {
+        std::shared_ptr<DGVertex> child = arc->dest();
+        // the vertex and the child must not be pinned to different locations
+        const bool location_conflict =
+            ((vptr)->address_set() || (vptr)->symbol_set()) &&
+            (child->address_set() || child->symbol_set());
+        if (!location_conflict && !child->is_a_target() &&
+            child->num_entry_arcs() == 1) {
+          if ((vptr)->address_set()) child->set_address((vptr)->address());
+          vptr->refer_this_to(child);
+          vptr->reset_symbol();
+          continue;
+        }
+      }
+    }
+
     // Is the exit arc DGArcDirect?
     {
       std::shared_ptr<DGArcDirect> arc_cast =
@@ -2270,6 +2296,9 @@ void DirectedGraph::update_func_names() {
   typedef vertices::iterator iter;
   for (iter v = stack_.begin(); v != stack_.end(); ++v) {
     const ver_ptr& vptr = vertex_ptr(*v);
+    // no code is generated for a vertex that refers to another (see
+    // handle_trivial_nodes()), hence no function is called for it
+    if ((vptr)->refers_to_another()) continue;
     // for every vertex with children
     if ((vptr)->num_exit_arcs() > 0) {
       // if it must be computed using a RR
diff --git a/src/bin/libint/extract.cc b/src/bin/libint/extract.cc
index 24b86252f..c395a5a1f 100644
--- a/src/bin/libint/extract.cc
+++ b/src/bin/libint/extract.cc
@@ -79,6 +79,10 @@ const ExtractExternSymbols::Symbols& ExtractExternSymbols::symbols() const {
 ////
 
 void ExtractRR::operator()(const VertexPtr& v) {
+  // no code is generated for a vertex that refers to another (e.g. a vertex
+  // whose RR is an alias, see DirectedGraph::handle_trivial_nodes()), hence
+  // neither is its RR used
+  if (v->refers_to_another()) return;
   if (v->num_exit_arcs() != 0) {
     std::shared_ptr<DGArc> arc = *(v->first_exit_arc());
     std::shared_ptr<DGArcRR> arc_rr =
diff --git a/src/bin/libint/master_ints_list.h b/src/bin/libint/master_ints_list.h
index 8ed891e8f..af528d1d3 100644
--- a/src/bin/libint/master_ints_list.h
+++ b/src/bin/libint/master_ints_list.h
@@ -159,16 +159,17 @@ typedef boost::mpl::list<
     σpeμσp_1_1_int, CMultipole_1_1_sh, CMultipole_1_1_int, CMultipole_1_1_sh_x,
     CMultipole_1_1_sh_y, CMultipole_1_1_sh_z, CMultipole_1_1_int_x,
     CMultipole_1_1_int_y, CMultipole_1_1_int_z, SMultipole_1_1_sh,
+    SMultipole_1_1_int
 #endif
-    SMultipole_1_1_int>
+    >
     MasterIntegralTypeList_1body_part;
 typedef boost::mpl::list<
     TwoPRep_11_11_sq, TwoPRep_11_11_int, Coulombσpσp_11_11_sq,
     Coulombσpσp_11_11_int, σpσpCoulombσpσp_11_11_sq, σpσpCoulombσpσp_11_11_int,
     σpCoulombσp_11_11_sq, σpCoulombσp_11_11_int, Coulombσp_11_11_sq,
-    Coulombσp_11_11_int, R12kG12_11_11_sq,
-    R12kG12_11_11_int, R12kR12lG12_11_11_sq, R12kR12lG12_11_11_int,
-    TiG12_11_11_sq, TiG12_11_11_int, G12TiG12_11_11_sq, G12TiG12_11_11_int,
+    Coulombσp_11_11_int, R12kG12_11_11_sq, R12kG12_11_11_int,
+    R12kR12lG12_11_11_sq, R12kR12lG12_11_11_int, TiG12_11_11_sq,
+    TiG12_11_11_int, G12TiG12_11_11_sq, G12TiG12_11_11_int,
     DivG12prime_xTx_11_11_sq, DivG12prime_xTx_11_11_int,
     DummySymmIntegral_11_11_sq, DummySymmIntegral_11_11_int>
     MasterIntegralTypeList_2body_part;
diff --git a/src/bin/libint/oper.h b/src/bin/libint/oper.h
index aec9a8867..100f5173f 100644
--- a/src/bin/libint/oper.h
+++ b/src/bin/libint/oper.h
@@ -30,6 +30,7 @@
 #include <vectorn.h>
 
 #include <boost/preprocessor/list/for_each.hpp>
+#include <stdexcept>
 #include <string>
 
 namespace libint2 {
@@ -304,7 +305,8 @@ struct σpVσp_Descr : public Contractable<σpVσp_Descr> {
 
   σpVσp_Descr() : quaternion_index_(0) {}
   σpVσp_Descr(int quaternion_index) : quaternion_index_(quaternion_index) {
-    assert(quaternion_index <= 3);
+    if (quaternion_index < 0 || quaternion_index > 3)
+      throw std::invalid_argument("quaternion_index out of range");
   }
 
   static const unsigned int max_key = 4;
@@ -345,7 +347,8 @@ struct σpeμσp_Descr : public Contractable<σpeμσp_Descr> {
 
   σpeμσp_Descr() : composite_index_(0) {}
   σpeμσp_Descr(int composite_index) : composite_index_(composite_index) {
-    assert(composite_index >= 0 && composite_index < 12);
+    if (composite_index < 0 || composite_index > 11)
+      throw std::invalid_argument("composite_index out of range");
   }
 
   static const unsigned int max_key = 12;
@@ -472,7 +475,8 @@ struct Coulombσpσp_Descr : public Contractable<Coulombσpσp_Descr> {
   Coulombσpσp_Descr() : quaternion_index_(0) {}
   Coulombσpσp_Descr(int quaternion_index)
       : quaternion_index_(quaternion_index) {
-    assert(quaternion_index <= 3);
+    if (quaternion_index < 0 || quaternion_index > 3)
+      throw std::invalid_argument("quaternion_index out of range");
   }
 
   static const unsigned int max_key = 4;
@@ -528,7 +532,8 @@ struct Coulombσp_Descr : public Contractable<Coulombσp_Descr> {
 
   Coulombσp_Descr() : cartesian_index_(0) {}
   Coulombσp_Descr(int cartesian_index) : cartesian_index_(cartesian_index) {
-    assert(cartesian_index >= 0 && cartesian_index <= 2);
+    if (cartesian_index < 0 || cartesian_index > 2)
+      throw std::invalid_argument("cartesian_index out of range");
   }
 
   static const unsigned int max_key = 3;
@@ -572,7 +577,8 @@ struct σpσpCoulombσpσp_Descr : public Contractable<σpσpCoulombσpσp_Descr
   σpσpCoulombσpσp_Descr() : component_index_(0) {}
   σpσpCoulombσpσp_Descr(int component_index)
       : component_index_(component_index) {
-    assert(component_index >= 0 && component_index <= 15);
+    if (component_index < 0 || component_index > 15)
+      throw std::invalid_argument("component_index out of range");
   }
 
   /// 16 components from tensor product of two independent spin spaces:
@@ -591,7 +597,7 @@ struct σpσpCoulombσpσp_Descr : public Contractable<σpσpCoulombσpσp_Descr
     };
     // clang-format on
     const auto qi = component_index();
-    if (qi > 15) abort();
+    if (qi < 0 || qi > 15) abort();
     return std::string("opop_coulomb_opop[") + labels[qi] + "]";
   }
   std::string label() const { return description(); }
@@ -642,7 +648,8 @@ struct σpCoulombσp_Descr : public Contractable<σpCoulombσp_Descr> {
 
   σpCoulombσp_Descr() : component_index_(0) {}
   σpCoulombσp_Descr(int component_index) : component_index_(component_index) {
-    assert(component_index >= 0 && component_index <= 8);
+    if (component_index < 0 || component_index > 8)
+      throw std::invalid_argument("component_index out of range");
   }
 
   static const unsigned int max_key = 9;
@@ -654,7 +661,7 @@ struct σpCoulombσp_Descr : public Contractable<σpCoulombσp_Descr> {
         "symtl_off_xy", "symtl_off_xz", "symtl_off_yz",
     };
     const auto ci = component_index();
-    if (ci > 8) abort();
+    if (ci < 0 || ci > 8) abort();
     return std::string("op_coulomb_op[") + labels[ci] + "]";
   }
   std::string label() const { return description(); }
diff --git a/src/bin/libint/rkb_fold_codegen.h b/src/bin/libint/rkb_fold_codegen.h
index 43360a91a..28a3b5d7d 100644
--- a/src/bin/libint/rkb_fold_codegen.h
+++ b/src/bin/libint/rkb_fold_codegen.h
@@ -107,9 +107,9 @@ void rkb_fold_make_children(Factory& factory, const F& a, const F& b,
   using libint2::prefactor::Scalar;
 
   const typename Factory::AuxIndexType zero_m(0u);
-  std::shared_ptr<DGVertex> first_child;
-  for (std::size_t t = 0; t != terms.size(); ++t) {
-    const RKBFoldTerm& term = terms[t];
+  std::vector<std::shared_ptr<DGVertex>> children;
+  children.reserve(terms.size());
+  for (const RKBFoldTerm& term : terms) {
     F fa{a};
     F fb{b};
     F fc{c};
@@ -118,28 +118,42 @@ void rkb_fold_make_children(Factory& factory, const F& a, const F& b,
     if (term.deriv[1] >= 0) fb.deriv().inc(term.deriv[1]);
     if (term.deriv[2] >= 0) fc.deriv().inc(term.deriv[2]);
     if (term.deriv[3] >= 0) fd.deriv().inc(term.deriv[3]);
-    const std::shared_ptr<DGVertex> child =
-        factory.make_child(fa, fb, fc, fd, zero_m);
-    if (!simple) continue;
+    children.push_back(factory.make_child(fa, fb, fc, fd, zero_m));
+  }
+  if (!simple) return;
 
-    if (t == 0) {
-      first_child = child;
-      // a bare child is not an expression: a leading +child is combined with
-      // the second term below, any other leading term is Scalar(coef) * child
-      if (term.coef != 1 || terms.size() == 1) expr = Scalar(term.coef) * child;
-    } else {
-      const std::shared_ptr<DGVertex> lhs =
-          (t == 1 && terms[0].coef == 1) ? first_child
-                                         : std::shared_ptr<DGVertex>(expr);
-      if (term.coef == 1)
-        expr = lhs + child;
-      else if (term.coef == -1)
-        expr = lhs - child;
-      else
-        expr = lhs + Scalar(term.coef) * child;
+  // the order of the children (= the order of the terms) is fixed by
+  // generate_code(), but the expression may sum them in any order: lead with
+  // a +1 term, if any, so that no Scalar(-1)*child is needed (a bare child is
+  // not an expression, hence a lone +1 term is Scalar(1)*child, which is free)
+  std::size_t lead = 0;
+  for (std::size_t t = 0; t != terms.size(); ++t)
+    if (terms[t].coef == 1) {
+      lead = t;
+      break;
+    }
+  std::shared_ptr<DGVertex> lhs = children[lead];
+  if (terms[lead].coef != 1) {
+    expr = Scalar(terms[lead].coef) * children[lead];
+    lhs = expr;
+    nflops += 1;
+  } else if (terms.size() == 1) {
+    expr = Scalar(1) * children[lead];
+    return;
+  }
+  for (std::size_t t = 0; t != terms.size(); ++t) {
+    if (t == lead) continue;
+    const auto coef = terms[t].coef;
+    if (coef == 1)
+      expr = lhs + children[t];
+    else if (coef == -1)
+      expr = lhs - children[t];
+    else {
+      expr = lhs + Scalar(coef) * children[t];
       nflops += 1;
     }
-    if (std::abs(term.coef) != 1) nflops += 1;
+    nflops += 1;
+    lhs = expr;
   }
 }
 
diff --git a/src/bin/libint/rr.h b/src/bin/libint/rr.h
index 1f46b65c1..bacc8172f 100644
--- a/src/bin/libint/rr.h
+++ b/src/bin/libint/rr.h
@@ -125,6 +125,13 @@ class RecurrenceRelation
      inlined and optimized.
   */
   virtual bool is_simple() const = 0;
+  /**
+     Returns true if the target is identical, element by element and in the
+     same layout, to the only child; then the vertex of the target is made to
+     refer to the child (see DirectedGraph::handle_trivial_nodes()), hence no
+     code is generated for this recurrence relation.
+  */
+  virtual bool is_alias() const { return false; }
   /**
      Returns true is the type of target and all children are exactly the same
   */
diff --git a/src/bin/libint/test.cc b/src/bin/libint/test.cc
index 9f9d0c061..ca47e3d0d 100644
--- a/src/bin/libint/test.cc
+++ b/src/bin/libint/test.cc
@@ -53,6 +53,7 @@ void test5();
 void test6();
 void test7();
 void test8();
+void test9();
 void test_cgshell_iter(const CGShell& sh);
 
 template <class Callback>
@@ -88,6 +89,7 @@ int main(int argc, char* argv[]) {
          << "  WARNING! Caught a standard exception:" << endl
          << "    " << a.what() << endl
          << endl;
+    return 1;
   }
   return 0;
 }
@@ -148,6 +150,9 @@ int try_main(int argc, char* argv[]) {
 #if 1
   RunTest(test8, "contracted G12 integral build");
 #endif
+#if 1
+  RunTest(test9, "DerivGaussV2 relation");
+#endif
 
   return 0;
 }
@@ -517,4 +522,148 @@ void test8() {
   cparams->contracted_targets(contracted_targets_old_value);
 }
 
+// DerivGaussV2: the children of direct differentiation and of translational
+// invariance, and the labels by which the RRStack shares the generated code
+void test9() {
+  typedef TwoPRep_11_11_sq IType;
+  auto check = [](bool ok, const std::string& what) {
+    if (!ok) throw ProgrammingError("DerivGaussV2 test: " + what);
+  };
+  // @return whether child @p i of @p rr is @p iset
+  auto child_is = [](const std::shared_ptr<RecurrenceRelation>& rr,
+                     unsigned int i, const std::shared_ptr<IType>& iset) {
+    return rr->rr_child(i) == std::static_pointer_cast<DGVertex, IType>(iset);
+  };
+  // @return shell of angular momentum @p l (primitive, unless @p contracted)
+  // with derivative quanta @p dx, @p dy, @p dz
+  auto shell = [](unsigned int l, unsigned int dx = 0, unsigned int dy = 0,
+                  unsigned int dz = 0, bool contracted = false) {
+    CGShell sh(l);
+    if (contracted)
+      sh.contract();
+    else
+      sh.uncontract();
+    if (dx) sh.deriv().inc(0, dx);
+    if (dy) sh.deriv().inc(1, dy);
+    if (dz) sh.deriv().inc(2, dz);
+    return sh;
+  };
+  const mType m0(0u);
+
+  // direct differentiation of a (particle 0, bra), in x:
+  //   ∂_Ax (p_x s|s s) = 2α (d s|s s) - 1 (s s|s s)
+  {
+    auto target =
+        IType::Instance(shell(1, 1), shell(0), shell(0), shell(0), m0);
+    auto rr = DerivV2_a_11_TwoPRep_11_sh::Instance(target, 0);
+    check(rr != nullptr, "no relation for (p_x s|s s)");
+    check(rr->num_children() == 2, "(p_x s|s s) must have 2 children");
+    check(child_is(rr, 0,
+                   IType::Instance(shell(2), shell(0), shell(0), shell(0), m0)),
+          "child 0 of (p_x s|s s) must be (d s|s s)");
+    check(child_is(rr, 1,
+                   IType::Instance(shell(0), shell(0), shell(0), shell(0), m0)),
+          "child 1 of (p_x s|s s) must be (s s|s s)");
+    // s has no a-1 term
+    auto target_s =
+        IType::Instance(shell(0, 1), shell(0), shell(0), shell(0), m0);
+    auto rr_s = DerivV2_a_11_TwoPRep_11_sh::Instance(target_s, 0);
+    check(rr_s != nullptr && rr_s->num_children() == 1,
+          "(s_x s|s s) must have 1 child");
+    check(child_is(rr_s, 0,
+                   IType::Instance(shell(1), shell(0), shell(0), shell(0), m0)),
+          "child 0 of (s_x s|s s) must be (p s|s s)");
+    // no derivative quanta in y, no relation
+    check(DerivV2_a_11_TwoPRep_11_sh::Instance(target, 1) == nullptr,
+          "(p_x s|s s) must not be differentiated in y");
+    // direct differentiation applies to primitives only
+    auto target_c =
+        IType::Instance(shell(1, 1, 0, 0, true), shell(0, 0, 0, 0, true),
+                        shell(0, 0, 0, 0, true), shell(0, 0, 0, 0, true), m0);
+    check(DerivV2_a_11_TwoPRep_11_sh::Instance(target_c, 0) == nullptr,
+          "contracted (p_x s|s s) must not be differentiated directly");
+  }
+
+  // code sharing: the label depends on the differentiated shell and the
+  // direction only, not on the spectators
+  {
+    auto rr1 = DerivV2_a_11_TwoPRep_11_sh::Instance(
+        IType::Instance(shell(1, 1), shell(0), shell(0), shell(0), m0), 0);
+    auto rr2 = DerivV2_a_11_TwoPRep_11_sh::Instance(
+        IType::Instance(shell(1, 1), shell(2), shell(1), shell(1), m0), 0);
+    auto rr3 = DerivV2_a_11_TwoPRep_11_sh::Instance(
+        IType::Instance(shell(2, 1), shell(0), shell(0), shell(0), m0), 0);
+    check(rr1 != nullptr && rr2 != nullptr && rr3 != nullptr,
+          "no relation for (p_x s|s s), (p_x d|p p) or (d_x s|s s)");
+    check(rr1->label() == rr2->label(),
+          "(p_x s|s s) and (p_x d|p p) must share the code: " + rr1->label() +
+              " vs " + rr2->label());
+    check(rr1->label() != rr3->label(),
+          "(p_x s|s s) and (d_x s|s s) must not share the code");
+    // a shell with derivative quanta in x and y is differentiated in either
+    // direction, with different children: the labels must differ
+    auto target_xy =
+        IType::Instance(shell(1, 1, 1), shell(0), shell(0), shell(0), m0);
+    auto rr_x = DerivV2_a_11_TwoPRep_11_sh::Instance(target_xy, 0);
+    auto rr_y = DerivV2_a_11_TwoPRep_11_sh::Instance(target_xy, 1);
+    check(rr_x != nullptr && rr_y != nullptr,
+          "(p_xy s|s s) must be differentiated in x and y");
+    check(rr_x->label() != rr_y->label(),
+          "the x and y derivatives of (p_xy s|s s) must not share the code");
+  }
+
+  // translational invariance (particle 1, bra = c), also for contracted
+  // shells: ∂_Cx (a b|c d) = -∂_Ax (a b|c d) - ∂_Bx (a b|c d) - ∂_Dx (a b|c d)
+  {
+    typedef DerivGaussV2<IType, 1, InBra, 1, InBra> TIType;
+    const bool c = true;
+    auto sh = [&](unsigned int l, unsigned int dx = 0) {
+      return shell(l, dx, 0, 0, c);
+    };
+    auto target = IType::Instance(sh(1), sh(0), sh(1, 1), sh(1), m0);
+    auto rr = TIType::Instance(target, 0);
+    check(rr != nullptr, "no translational invariance relation");
+    check(rr->num_children() == 3,
+          "translational invariance must have 3 children");
+    check(child_is(rr, 0, IType::Instance(sh(1, 1), sh(0), sh(1), sh(1), m0)),
+          "child 0 must be (p_x s|p p)");
+    check(child_is(rr, 1, IType::Instance(sh(1), sh(0, 1), sh(1), sh(1), m0)),
+          "child 1 must be (p s_x|p p)");
+    check(child_is(rr, 2, IType::Instance(sh(1), sh(0), sh(1), sh(1, 1), m0)),
+          "child 2 must be (p s|p p_x)");
+    check(rr->label().rfind("DerivGaussV2TI", 0) == 0,
+          "translational invariance label must start with DerivGaussV2TI: " +
+              rr->label());
+  }
+
+  // build contracted shell sets, without unrolling, hence with calls to the
+  // generated DerivGaussV2 functions: differentiation of a (direct, shared by
+  // the spectators) and of c (translational invariance)
+  {
+    CGShell::set_contracted_default_value(true);
+    const bool contracted_targets_old_value = cparams->contracted_targets();
+    cparams->contracted_targets(true);
+    const bool c = true;
+    RunBuildTest<IType>(shell(1, 1, 0, 0, c), shell(1, 0, 0, 0, c),
+                        shell(0, 0, 0, 0, c), shell(1, 0, 0, 0, c), 0,
+                        use_quartets);
+    RunBuildTest<IType>(shell(1, 0, 0, 0, c), shell(0, 0, 0, 0, c),
+                        shell(1, 0, 1, 0, c), shell(0, 0, 0, 0, c), 0,
+                        use_quartets);
+    cparams->contracted_targets(contracted_targets_old_value);
+  }
+
+  // CR_11_Coulombσp_11 is an alias (target = child, see
+  // DirectedGraph::handle_trivial_nodes()), CR_11_Coulombσpσp_11 is not
+  {
+    for (int k = 0; k != 3; ++k) {
+      auto target = Coulombσp_11_11_sq::Instance(
+          shell(1), shell(0), shell(1), shell(1), m0, Coulombσp_Descr(k));
+      auto rr = CR_11_Coulombσp_11_sh::Instance(target, 0);
+      check(rr != nullptr && rr->num_children() == 1 && rr->is_alias(),
+            "CR_11_Coulombσp_11 must be an alias");
+    }
+  }
+}
+
 };  // namespace
diff --git a/src/lib/libint/features.cmake.in b/src/lib/libint/features.cmake.in
index 2d310e873..0d3fff762 100644
--- a/src/lib/libint/features.cmake.in
+++ b/src/lib/libint/features.cmake.in
@@ -3,6 +3,7 @@ LIBINT_ONEBODY_DERIV = @LIBINT_ONEBODY_DERIV@
 LIBINT_SUPPORTS_ERI = @LIBINT_SUPPORTS_ERI@
 LIBINT_ERI_MAX_AM = @Libint2_MAX_AM_ERI@
 LIBINT_ERI_DERIV = @LIBINT_ERI_DERIV@
+LIBINT_ERI_MAX_AM_DERIVS = [@Libint2_MAX_AM_ERI_DERIVS@]
 LIBINT_CONTRACTED_INTS = @LIBINT_CONTRACTED_INTS@
 LIBINT_SHELL_SET = @LIBINT_SHELL_SET@
 LIBINT_SHGSHELL_ORDERING = @LIBINT_SHGSHELL_ORDERING@