From 996585cb6a0921b01901f1243d5577e8adf81f27 Mon Sep 17 00:00:00 2001
From: Jonathon Misiewicz <jonathon.paul.misiewicz@emory.edu>
Date: Mon, 27 Nov 2023 15:31:11 -0500
Subject: [PATCH] Manage the Boys number.

---
 include/libint2/engine.h      |  5 +++++
 include/libint2/engine.impl.h | 27 +++++++++++++++++++++++----
 2 files changed, 28 insertions(+), 4 deletions(-)

diff --git a/include/libint2/engine.h b/include/libint2/engine.h
index 8d99e519a..4cf64c3ff 100644
--- a/include/libint2/engine.h
+++ b/include/libint2/engine.h
@@ -187,6 +187,9 @@ struct default_operator_traits {
   } oper_params_type;
   static oper_params_type default_params() { return oper_params_type{}; }
   static constexpr auto nopers = 1u;
+  // N.B.: Below field means we *should* template specialize operator_traits for
+  // Operator::kinetic, but L2 doesn't use that anywhere.
+  static constexpr auto intrinsic_nderiv = 0u;
   struct _core_eval_type {
     template <typename... params>
     static std::shared_ptr<const _core_eval_type> instance(params...) {
@@ -222,6 +225,7 @@ template <>
 struct operator_traits<Operator::σpVσp>
     : public operator_traits<Operator::nuclear> {
   static constexpr auto nopers = 4;
+  static constexpr auto intrinsic_nderiv = 2;
 };
 
 template <>
@@ -1000,6 +1004,7 @@ class Engine {
   //-------
   unsigned int nparams() const;
   unsigned int nopers() const;
+  unsigned int intrinsic_nderiv() const;
   /// if Params == operator_traits<oper>::oper_params_type, will return
   /// \c any(params)
   /// else will set return \c any initialized with default value for
diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h
index 0f4fc8f13..119043a98 100644
--- a/include/libint2/engine.impl.h
+++ b/include/libint2/engine.impl.h
@@ -775,6 +775,19 @@ __libint2_engine_inline unsigned int Engine::nopers() const {
   assert(false && "missing case in switch");  // unreachable
   abort();
 }
+__libint2_engine_inline unsigned int Engine::intrinsic_nderiv() const {
+  switch (static_cast<int>(oper_)) {
+#define BOOST_PP_NBODYENGINE_MCR9(r, data, i, elem) \
+  case i:                                           \
+    return operator_traits<static_cast<Operator>(i)>::intrinsic_nderiv;
+    BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR9, _,
+                             BOOST_PP_NBODY_OPERATOR_LIST)
+    default:
+      break;
+  }
+  assert(false && "missing case in switch");  // unreachable
+  abort();
+}
 
 template <>
 __libint2_engine_inline any Engine::enforce_params_type<any>(
@@ -828,6 +841,10 @@ __libint2_engine_inline any Engine::enforce_params_type(
   return result;
 }
 
+/// @param[in] oper the Operator for which to return core-evaluator objects
+/// @return a core-integral evaluator and its scratch for @c oper
+/// @throw @c oper not registered in BOOST_PP_NBODY_OPERATOR_LIST
+/// @throw @c oper has no core_eval_type specialization in operator_traits.
 __libint2_engine_inline any Engine::make_core_eval_pack(Operator oper) const {
   any result;
   switch (static_cast<int>(oper)) {
@@ -835,11 +852,13 @@ __libint2_engine_inline any Engine::make_core_eval_pack(Operator oper) const {
   case i:                                                                    \
     result = libint2::detail::make_compressed_pair(                          \
         operator_traits<static_cast<Operator>(i)>::core_eval_type::instance( \
-            braket_rank() * lmax_ + deriv_order_,                            \
+            braket_rank() * lmax_ + deriv_order_ +                           \
+            operator_traits<static_cast<Operator>(i)>::intrinsic_nderiv,     \
             std::numeric_limits<scalar_type>::epsilon()),                    \
         libint2::detail::CoreEvalScratch<                                    \
             operator_traits<static_cast<Operator>(i)>::core_eval_type>(      \
-            braket_rank() * lmax_ + deriv_order_));                          \
+            braket_rank() * lmax_ + deriv_order_ +                           \
+            operator_traits<static_cast<Operator>(i)>::intrinsic_nderiv));   \
     assert(any_cast<detail::core_eval_pack_type<static_cast<Operator>(i)>>(  \
                &result) != nullptr);                                         \
     break;
@@ -1057,9 +1076,9 @@ __libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata, const
                      primdata.PC_y[0] * primdata.PC_y[0] +
                      primdata.PC_z[0] * primdata.PC_z[0];
     const scalar_type U = gammap * PC2;
-    const auto mmax = s1.contr[0].l + s2.contr[0].l + deriv_order_;
+    const auto mmax = s1.contr[0].l + s2.contr[0].l + deriv_order_ + intrinsic_nderiv();
     auto* fm_ptr = &(primdata.LIBINT_T_S_ELECPOT_S(0)[0]);
-    if (oper_ == Operator::nuclear) {
+    if (oper_ == Operator::nuclear || oper_ == Operator::σpVσp) {
       const auto& fm_engine_ptr =
           any_cast<const detail::core_eval_pack_type<Operator::nuclear>&>(core_eval_pack_)
           .first();