engine.h

/*
 *  Copyright (C) 2004-2021 Edward F. Valeev
 *
 *  This file is part of Libint.
 *
 *  Libint is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  Libint is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with Libint.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef _libint2_src_lib_libint_engine_h_
#define _libint2_src_lib_libint_engine_h_
 
#ifndef LIBINT2_DOES_NOT_INLINE_ENGINE
# define __libint2_engine_inline inline
#else
# define __libint2_engine_inline
#endif
 
#include <libint2/util/cxxstd.h>
#if LIBINT2_CPLUSPLUS_STD < 2011
# error "libint2/engine.h requires C++11 support"
#endif
 
#include <algorithm>
#include <array>
#include <cstring>
#include <cstdio>
#include <functional>
#include <iostream>
#include <limits>
#include <map>
#include <memory>
#include <tuple>
#include <utility>
#include <vector>
 
#include <libint2/cxxapi.h>
#include <libint2/boys_fwd.h>
#include <libint2/shell.h>
#include <libint2/solidharmonics.h>
#include <libint2/cartesian.h>
#include <libint2/util/any.h>
#include <libint2/util/array_adaptor.h>
#include <libint2/util/intpart_iter.h>
#include <libint2/util/compressed_pair.h>
#include <libint2/util/timer.h>
#include <libint2/cgshellinfo.h>
#include <libint2/braket.h>
 
// the engine will be profiled by default if library was configured with
// --enable-profile
#ifdef LIBINT2_PROFILE
#define LIBINT2_ENGINE_TIMERS
// uncomment if want to profile each integral class
#define LIBINT2_ENGINE_PROFILE_CLASS
#endif
// uncomment if want to profile the engine even if library was configured
// without --enable-profile
//#  define LIBINT2_ENGINE_TIMERS
 
namespace libint2 {
 
typedef std::vector<std::pair<double, double>> ContractedGaussianGeminal;
 
constexpr size_t num_geometrical_derivatives(size_t ncenter,
                                             size_t deriv_order) {
  return (deriv_order > 0)
             ? (num_geometrical_derivatives(ncenter, deriv_order - 1) *
                (3 * ncenter + deriv_order - 1)) /
                   deriv_order
             : 1;
}
 
template <typename T, unsigned N>
__libint2_engine_inline typename std::remove_all_extents<T>::type* to_ptr1(T (&a)[N]);
 
 
enum class Operator {
  overlap = 0,
  kinetic,
  nuclear,
  erf_nuclear,
  erfc_nuclear,
  emultipole1,
  emultipole2,
  emultipole3,
  sphemultipole,
  delta,
  coulomb,
  r12_m1 = coulomb,
  cgtg,
  cgtg_x_coulomb,
  delcgtg2,
  r12,
  r12_1 = r12,
  erf_coulomb,
  erfc_coulomb,
  stg,
  stg_x_coulomb,
  yukawa = stg_x_coulomb,
  // do not modify this
  invalid = -1,
  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!keep this updated!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  first_1body_oper = overlap,
  last_1body_oper = sphemultipole,
  first_2body_oper = delta,
  last_2body_oper = stg_x_coulomb,
  first_oper = first_1body_oper,
  last_oper = last_2body_oper
};
 
inline constexpr int rank(Operator oper) {
  return (oper >= Operator::first_1body_oper &&
          oper <= Operator::last_1body_oper)
          ? 1 :
         ((oper >= Operator::first_2body_oper &&
           oper <= Operator::last_2body_oper)
           ? 2 : throw std::logic_error("rank(Operator): invalid operator given"));
}
 
namespace detail {
struct default_operator_traits {
  typedef struct {
  } oper_params_type;
  static oper_params_type default_params() { return oper_params_type{}; }
  static constexpr auto nopers = 1u;
  struct _core_eval_type {
    template <typename... params>
    static std::shared_ptr<const _core_eval_type> instance(params...) {
      return nullptr;
    }
  };
  using core_eval_type = const _core_eval_type;
};
}  // namespace detail
 
template <Operator Op>
struct operator_traits : public detail::default_operator_traits {};
 
template <>
struct operator_traits<Operator::nuclear>
    : public detail::default_operator_traits {
  typedef std::vector<std::pair<scalar_type, std::array<scalar_type, 3>>>
      oper_params_type;
  static oper_params_type default_params() { return oper_params_type{}; }
#ifndef LIBINT_USER_DEFINED_REAL
  typedef const libint2::FmEval_Chebyshev7<scalar_type> core_eval_type;
#else
  typedef const libint2::FmEval_Reference<scalar_type> core_eval_type;
#endif
};
 
template <>
struct operator_traits<Operator::erf_nuclear>
    : public detail::default_operator_traits {
  typedef std::tuple<
      scalar_type, typename operator_traits<Operator::nuclear>::oper_params_type>
      oper_params_type;
  static oper_params_type default_params() {
    return std::make_tuple(0,operator_traits<Operator::nuclear>::default_params());
  }
  typedef const libint2::GenericGmEval<libint2::os_core_ints::erf_coulomb_gm_eval<scalar_type>>
      core_eval_type;
};
 
template <>
struct operator_traits<Operator::erfc_nuclear>
    : public detail::default_operator_traits {
  typedef typename operator_traits<Operator::erf_nuclear>::oper_params_type oper_params_type;
  static oper_params_type default_params() {
    return std::make_tuple(0,operator_traits<Operator::nuclear>::default_params());
  }
  typedef const libint2::GenericGmEval<libint2::os_core_ints::erfc_coulomb_gm_eval<scalar_type>>
      core_eval_type;
};
 
template <>
struct operator_traits<Operator::emultipole1>
    : public detail::default_operator_traits {
  typedef std::array<double, 3> oper_params_type;
  static oper_params_type default_params() {
    return oper_params_type{{0,0,0}};
  }
  static constexpr auto nopers = 4u;  
};
template <>
struct operator_traits<Operator::emultipole2>
    : public operator_traits<Operator::emultipole1> {
  static constexpr auto nopers =
      operator_traits<Operator::emultipole1>::nopers +
      6;  
};
template <>
struct operator_traits<Operator::emultipole3>
    : public operator_traits<Operator::emultipole1> {
  static constexpr auto nopers =
      operator_traits<Operator::emultipole2>::nopers + 10;
};
template <>
struct operator_traits<Operator::sphemultipole>
    : public operator_traits<Operator::emultipole1> {
  static constexpr auto nopers =
      (MULTIPOLE_MAX_ORDER + 1) * (MULTIPOLE_MAX_ORDER + 1);
};
 
template <>
struct operator_traits<Operator::coulomb>
    : public detail::default_operator_traits {
#ifndef LIBINT_USER_DEFINED_REAL
  typedef const libint2::FmEval_Chebyshev7<scalar_type> core_eval_type;
#else
  typedef const libint2::FmEval_Reference<scalar_type> core_eval_type;
#endif
};
namespace detail {
template <int K>
struct cgtg_operator_traits : public detail::default_operator_traits {
  typedef libint2::GaussianGmEval<scalar_type, K> core_eval_type;
  typedef ContractedGaussianGeminal oper_params_type;
};
}  // namespace detail
template <>
struct operator_traits<Operator::cgtg>
    : public detail::cgtg_operator_traits<0> {};
template <>
struct operator_traits<Operator::cgtg_x_coulomb>
    : public detail::cgtg_operator_traits<-1> {};
template <>
struct operator_traits<Operator::delcgtg2>
    : public detail::cgtg_operator_traits<2> {};
 
template <>
struct operator_traits<Operator::delta>
    : public detail::default_operator_traits {
  typedef const libint2::GenericGmEval<libint2::os_core_ints::delta_gm_eval<scalar_type>>
      core_eval_type;
};
 
template <>
struct operator_traits<Operator::r12>
    : public detail::default_operator_traits {
  typedef const libint2::GenericGmEval<libint2::os_core_ints::r12_xx_K_gm_eval<scalar_type, 1>>
      core_eval_type;
};
 
template <>
struct operator_traits<Operator::erf_coulomb>
    : public detail::default_operator_traits {
  typedef scalar_type oper_params_type;
  static oper_params_type default_params() {
    return oper_params_type{0};
  }
  typedef const libint2::GenericGmEval<libint2::os_core_ints::erf_coulomb_gm_eval<scalar_type>>
      core_eval_type;
};
template <>
struct operator_traits<Operator::erfc_coulomb>
    : public detail::default_operator_traits {
  typedef scalar_type oper_params_type;
  static oper_params_type default_params() {
    return oper_params_type{0};
  }
  typedef const libint2::GenericGmEval<libint2::os_core_ints::erfc_coulomb_gm_eval<scalar_type>>
      core_eval_type;
};
 
template <>
struct operator_traits<Operator::stg>
    : public detail::default_operator_traits {
  typedef scalar_type oper_params_type;
  static oper_params_type default_params() {
    return oper_params_type{0};
  }
  typedef const libint2::TennoGmEval<scalar_type>
      core_eval_type;
};
 
template <>
struct operator_traits<Operator::stg_x_coulomb>
    : public detail::default_operator_traits {
  typedef scalar_type oper_params_type;
  static oper_params_type default_params() {
    return oper_params_type{0};
  }
  typedef const libint2::TennoGmEval<scalar_type>
      core_eval_type;
};
 
libint2::any
default_params(const Operator& oper);
 
namespace detail {
template <typename core_eval_type>
using __core_eval_pack_type =
    compressed_pair<std::shared_ptr<core_eval_type>,
                    libint2::detail::CoreEvalScratch<core_eval_type>>;
template <Operator Op>
using core_eval_pack_type =
    __core_eval_pack_type<typename operator_traits<Op>::core_eval_type>;
}
 
#define BOOST_PP_NBODY_BRAKET_MAX_INDEX 4
 
inline constexpr int rank(BraKet braket) {
  return (braket==BraKet::x_x || braket==BraKet::xs_xs) ? 2 :
         ((braket==BraKet::xs_xx || braket==BraKet::xx_xs) ? 3 :
          ((braket==BraKet::xx_xx) ? 4 : throw std::logic_error("rank(BraKet): invalid braket given")));
}
 
inline constexpr BraKet default_braket(Operator oper) {
  return (rank(oper)==1) ? BraKet::x_x :
         ((rank(oper)==2) ? BraKet::xx_xx : throw std::logic_error("default_braket(Operator): invalid operator given"));
}
 
constexpr size_t nopers_2body = static_cast<int>(Operator::last_2body_oper) -
                                static_cast<int>(Operator::first_2body_oper) +
                                1;
constexpr size_t nbrakets_2body = static_cast<int>(BraKet::last_2body_braket) -
                                  static_cast<int>(BraKet::first_2body_braket) +
                                  1;
constexpr size_t nderivorders_2body = LIBINT2_MAX_DERIV_ORDER + 1;
 
class Engine {
 private:
  typedef struct {
  } empty_pod;
 
 public:
 
  static constexpr auto max_ntargets =
      std::extent<decltype(std::declval<Libint_t>().targets), 0>::value;
  using target_ptr_vec =
      std::vector<const value_type*, detail::ext_stack_allocator<const value_type*, max_ntargets>>;
 
  Engine()
      : oper_(Operator::invalid),
        braket_(BraKet::invalid),
        primdata_(),
        stack_size_(0),
        lmax_(-1),
        deriv_order_(0),
        cartesian_shell_normalization_(CartesianShellNormalization::standard),
        scale_(1) {
    set_precision(std::numeric_limits<scalar_type>::epsilon());
  }
 
  // clang-format off
 
  // clang-format on
  template <typename Params = empty_pod>
  Engine(Operator oper, size_t max_nprim, int max_l, int deriv_order = 0,
         scalar_type precision = std::numeric_limits<scalar_type>::epsilon(),
         Params params = empty_pod(), BraKet braket = BraKet::invalid,
         ScreeningMethod screening_method = default_screening_method())
      : oper_(oper),
        braket_(braket),
        primdata_(),
        spbra_(max_nprim),
        spket_(max_nprim),
        stack_size_(0),
        lmax_(max_l),
        deriv_order_(deriv_order),
        screening_method_(screening_method),
        cartesian_shell_normalization_(CartesianShellNormalization::standard),
        scale_(1),
        params_(enforce_params_type(oper, params)) {
    set_precision(precision);
    assert(max_nprim > 0);
    initialize(max_nprim);
    core_eval_pack_ = make_core_eval_pack(oper);  // must follow initialize() to
                                                  // ensure default braket_ has
                                                  // been set
    init_core_ints_params(params_);
  }
 
  // intel does not accept "move ctor = default"
  Engine(Engine&& other)
      : oper_(other.oper_),
        braket_(other.braket_),
        primdata_(std::move(other.primdata_)),
        spbra_(std::move(other.spbra_)),
        spket_(std::move(other.spket_)),
        stack_size_(other.stack_size_),
        lmax_(other.lmax_),
        hard_lmax_(other.hard_lmax_),
        hard_default_lmax_(other.hard_default_lmax_),
        deriv_order_(other.deriv_order_),
        precision_(other.precision_),
        ln_precision_(other.ln_precision_),
        screening_method_(other.screening_method_),
        cartesian_shell_normalization_(other.cartesian_shell_normalization_),
        scale_(other.scale_),
        core_eval_pack_(std::move(other.core_eval_pack_)),
        params_(std::move(other.params_)),
        core_ints_params_(std::move(other.core_ints_params_)),
        targets_(std::move(other.targets_)),
        set_targets_(other.set_targets_),
        scratch_(std::move(other.scratch_)),
        scratch2_(other.scratch2_),
        buildfnptrs_(other.buildfnptrs_) {
    // leave other in an unusable (but valid) state
    other.oper_ = Operator::invalid;
    other.braket_ = BraKet::invalid;
    other.scratch2_ = nullptr;
  }
 
  Engine(const Engine& other)
      : oper_(other.oper_),
        braket_(other.braket_),
        primdata_(other.primdata_.size()),
        spbra_(other.spbra_),
        spket_(other.spket_),
        stack_size_(other.stack_size_),
        lmax_(other.lmax_),
        deriv_order_(other.deriv_order_),
        precision_(other.precision_),
        ln_precision_(other.ln_precision_),
        screening_method_(other.screening_method_),
        cartesian_shell_normalization_(other.cartesian_shell_normalization_),
        scale_(other.scale_),
        core_eval_pack_(other.core_eval_pack_),
        params_(other.params_),
        core_ints_params_(other.core_ints_params_) {
    initialize();
  }
 
  ~Engine() { finalize(); }
 
  Engine& operator=(Engine&& other) {
    oper_ = other.oper_;
    braket_ = other.braket_;
    primdata_ = std::move(other.primdata_);
    spbra_ = std::move(other.spbra_);
    spket_ = std::move(other.spket_);
    stack_size_ = other.stack_size_;
    lmax_ = other.lmax_;
    hard_lmax_ = other.hard_lmax_;
    hard_default_lmax_ = other.hard_default_lmax_;
    deriv_order_ = other.deriv_order_;
    precision_ = other.precision_;
    ln_precision_ = other.ln_precision_;
    screening_method_ = other.screening_method_;
    cartesian_shell_normalization_ = other.cartesian_shell_normalization_;
    scale_ = other.scale_;
    core_eval_pack_ = std::move(other.core_eval_pack_);
    params_ = std::move(other.params_);
    core_ints_params_ = std::move(other.core_ints_params_);
    targets_ = std::move(other.targets_);
    set_targets_ = other.set_targets_;
    scratch_ = std::move(other.scratch_);
    scratch2_ = other.scratch2_;
    buildfnptrs_ = other.buildfnptrs_;
    // leave other in an unusable state
    other.oper_ = Operator::invalid;
    other.braket_ = BraKet::invalid;
    other.scratch2_ = nullptr;
    return *this;
  }
 
  Engine& operator=(const Engine& other) {
    oper_ = other.oper_;
    braket_ = other.braket_;
    primdata_.resize(other.primdata_.size());
    spbra_ = other.spbra_;
    spket_ = other.spket_;
    stack_size_ = other.stack_size_;
    lmax_ = other.lmax_;
    deriv_order_ = other.deriv_order_;
    precision_ = other.precision_;
    ln_precision_ = other.ln_precision_;
    screening_method_ = other.screening_method_;
    cartesian_shell_normalization_ = other.cartesian_shell_normalization_;
    scale_ = other.scale_;
    core_eval_pack_ = other.core_eval_pack_;
    params_ = other.params_;
    core_ints_params_ = other.core_ints_params_;
    initialize();
    return *this;
  }
 
  int operator_rank() const { return rank(oper_); }
 
  int braket_rank() const { return rank(braket_); }
 
  Operator oper() const { return oper_; }
 
  BraKet braket() const { return braket_; }
 
  int deriv_order() const { return deriv_order_; }
 
  Engine& set(Operator new_oper) {
    if (oper_ != new_oper) {
      oper_ = new_oper;
      braket_ = default_braket(oper_);
      params_ = default_params(oper_);
      initialize();
      core_eval_pack_ = make_core_eval_pack(oper_);  // must follow initialize() to
                                                     // ensure default braket_ has
                                                     // been set
    }
    return *this;
  }
 
  Engine& set(BraKet new_braket) {
    if (braket_ != new_braket) {
      braket_ = new_braket;
      initialize();
    }
    return *this;
  }
 
  template <typename Params>
  Engine& set_params(const Params& params) {
    params_ = params;
    init_core_ints_params(params_);
    reset_scratch();
    return *this;
  }
 
  std::size_t max_nprim() const {
    assert(spbra_.primpairs.size() == spket_.primpairs.size());
    return static_cast<std::size_t>(std::sqrt(spbra_.primpairs.size()));
  }
 
  Engine& set_max_nprim(std::size_t n) {
    if (n*n > spbra_.primpairs.size()) {
      spbra_.resize(n);
      spket_.resize(n);
      initialize(n);
    }
    return *this;
  }
 
  std::size_t max_l() const {
    return lmax_;
  }
 
  Engine& set_max_l(std::size_t L) {
    if (L >= static_cast<std::size_t>(lmax_)) {
      lmax_ = L;
      initialize();
    }
    return *this;
  }
 
  const target_ptr_vec& results() const { return targets_; }
 
  unsigned int nshellsets() const { return targets_.size(); }
 
 
  template <typename... ShellPack>
  __libint2_engine_inline const target_ptr_vec& compute(
      const libint2::Shell& first_shell, const ShellPack&... rest_of_shells);
 
  __libint2_engine_inline const target_ptr_vec& compute1(
      const libint2::Shell& s1, const libint2::Shell& s2);
 
  template <Operator oper, BraKet braket, size_t deriv_order>
  __libint2_engine_inline const target_ptr_vec& compute2(const Shell& bra1,
                                                         const Shell& bra2,
                                                         const Shell& ket1,
                                                         const Shell& ket2,
                                                         const ShellPair* spbra = nullptr,
                                                         const ShellPair* spket = nullptr);
 
  typedef const target_ptr_vec& (Engine::*compute2_ptr_type)(const Shell& bra1,
                                                             const Shell& bra2,
                                                             const Shell& ket1,
                                                             const Shell& ket2,
                                                             const ShellPair* spbra,
                                                             const ShellPair* spket);
 
  // clang-format off
  // clang-format on
  Engine& set_precision(scalar_type prec) {
    if (prec <= 0.) {
      precision_ = 0.;
      ln_precision_ = std::numeric_limits<scalar_type>::lowest();
      return *this;
    }
    if (prec > 0.) {  // split single if/else to avoid speculative execution of else branch on Apple M1 with apple clang 13.0.0
      precision_ = prec;
      using std::log;
      ln_precision_ = log(precision_);
      return *this;
    }
    abort();
  }
  scalar_type precision() const { return precision_; }
 
  Engine& set(ScreeningMethod screening_method) { screening_method_ = screening_method; return *this; }
 
  ScreeningMethod screening_method() const { return screening_method_; }
 
  CartesianShellNormalization cartesian_shell_normalization() const {
    return cartesian_shell_normalization_;
  }
 
  Engine& set(CartesianShellNormalization norm) {
    cartesian_shell_normalization_ = norm;
    return *this;
  }
 
  scalar_type prescaled_by() const {
    return scale_;
  }
 
  Engine& prescale_by(scalar_type sc) {
    scale_ = sc;
    return *this;
  }
 
  void print_timers(std::ostream& os = std::cout) {
#ifdef LIBINT2_ENGINE_TIMERS
    os << "timers: prereq = " << timers.read(0);
#ifndef LIBINT2_PROFILE  // if libint's profiling was on, engine's build timer
                         // will include its overhead
    // do not report it, detailed profiling from libint will be printed below
    os << " build = " << timers.read(1);
#endif
    os << " tform = " << timers.read(2) << std::endl;
#endif
#ifdef LIBINT2_PROFILE
    os << "build timers: hrr = " << primdata_[0].timers->read(0)
       << " vrr = " << primdata_[0].timers->read(1) << std::endl;
#endif
#ifdef LIBINT2_ENGINE_TIMERS
#ifdef LIBINT2_ENGINE_PROFILE_CLASS
    for (const auto& p : class_profiles) {
      char buf[1024];
      std::snprintf(buf, sizeof(buf), "{\"%s\", %10.5lf, %10.5lf, %10.5lf, %10.5lf, %ld, %ld},",
             p.first.to_string().c_str(), p.second.prereqs, p.second.build_vrr,
             p.second.build_hrr, p.second.tform, p.second.nshellset,
             p.second.nprimset);
      os << buf << std::endl;
    }
#endif
#endif
  }
 
  class lmax_exceeded : virtual public std::logic_error {
   public:
    lmax_exceeded(const char* task_name, size_t lmax_limit,
                  size_t lmax_requested)
        : std::logic_error(
              "Engine::lmax_exceeded -- angular momentum limit exceeded"),
          lmax_limit_(lmax_limit),
          lmax_requested_(lmax_requested) {
      strncpy(task_name_, task_name, 64);
      task_name_[64] = '\0';
    }
    ~lmax_exceeded() noexcept {}
 
    const char* task_name() const {
      return static_cast<const char*>(task_name_);
    }
    size_t lmax_limit() const { return lmax_limit_; }
    size_t lmax_requested() const { return lmax_requested_; }
 
   private:
    char task_name_[65];
    size_t lmax_limit_;
    size_t lmax_requested_;
  };
 
  class using_default_initialized : virtual public std::logic_error {
   public:
    using_default_initialized()
        : std::logic_error(
        "Engine::using_default_initialized -- attempt to use a default-initialized Engine") {
    }
    ~using_default_initialized() noexcept {}
  };
 
 private:
  Operator oper_;
  BraKet braket_;
  std::vector<Libint_t> primdata_;
  ShellPair spbra_, spket_;
  size_t stack_size_;  // amount allocated by libint2_init_xxx in
                       // primdata_[0].stack
  int lmax_;
  int hard_lmax_;  // max L supported by library for this operator type (i.e. LIBINT2_MAX_AM_<task><deriv_order>) + 1
  int hard_default_lmax_;  // max L supported by library for default operator type
                           // (i.e. LIBINT2_MAX_AM_default<deriv_order>) + 1
                           // this is only set and used if LIBINT2_CENTER_DEPENDENT_MAX_AM_<task><deriv_order> == 1
  int deriv_order_;
  scalar_type precision_;
  scalar_type ln_precision_;
  ScreeningMethod screening_method_ = ScreeningMethod::Invalid;
 
  // specifies the normalization convention for Cartesian Gaussians
  CartesianShellNormalization cartesian_shell_normalization_;
 
  // the target integrals are scaled by this factor (of course it is actually the core integrals that are scaled)
  scalar_type scale_;
 
  any core_eval_pack_;
 
  any params_;  // operator params
  any core_ints_params_;
  void init_core_ints_params(const any& params);
 
  target_ptr_vec targets_;
  bool set_targets_;
 
  std::vector<value_type>
      scratch_;       // for transposes and/or transforming to solid harmonics
  value_type* scratch2_;  // &scratch_[0] points to the first block large enough to
                      // hold all target ints
  // scratch2_ points to second such block. It could point into scratch_ or at
  // primdata_[0].stack
  typedef void (*buildfnptr_t)(const Libint_t*);
  buildfnptr_t* buildfnptrs_;
 
  unsigned int compute_nshellsets() const {
    const unsigned int num_operator_geometrical_derivatives =
        (oper_ == Operator::nuclear || oper_ == Operator::erf_nuclear ||
         oper_ == Operator::erfc_nuclear)
            ? this->nparams()
            : 0;
    const auto ncenters = braket_rank() + num_operator_geometrical_derivatives;
    return nopers() * num_geometrical_derivatives(ncenters, deriv_order_);
  }
 
  void reset_scratch() {
    const auto nshsets = compute_nshellsets();
    targets_.resize(nshsets);
    set_targets_ = (&targets_[0] != const_cast<const value_type**>(primdata_[0].targets));
    const auto ncart_max = (lmax_ + 1) * (lmax_ + 2) / 2;
    const auto target_shellset_size =
        nshsets * std::pow(ncart_max, braket_rank());
    // need to be able to hold 2 sets of target shellsets: the worst case occurs
    // when dealing with
    // 1-body Coulomb ints derivatives ... have 2+natom 1st-order derivative sets
    // (and many more of 2nd and higher) that are stored in scratch
    // then need to transform to solids. To avoid copying back and forth make
    // sure that there is enough
    // room to transform all ints and save them in correct order in single pass
    const auto need_extra_large_scratch = stack_size_ < target_shellset_size;
    scratch_.resize(need_extra_large_scratch ? 2 * target_shellset_size
                                             : target_shellset_size);
    scratch2_ = need_extra_large_scratch ? &scratch_[target_shellset_size]
                                         : primdata_[0].stack;
  }
 
  __libint2_engine_inline void compute_primdata(Libint_t& primdata,
                                                const Shell& s1,
                                                const Shell& s2, size_t p1,
                                                size_t p2, size_t oset);
 
public:
  __libint2_engine_inline const std::vector<Engine::compute2_ptr_type>&
  compute2_ptrs() const;
 
private:
  // max_nprim=0 avoids resizing primdata_
  __libint2_engine_inline void initialize(size_t max_nprim = 0);
  // generic _initializer
  __libint2_engine_inline void _initialize();
 
  void finalize() {
    if (primdata_.size() != 0) {
      libint2_cleanup_default(&primdata_[0]);
    }
  }  // finalize()
 
  //-------
  // utils
  //-------
  unsigned int nparams() const;
  unsigned int nopers() const;
  template <typename Params>
  static any enforce_params_type(
      Operator oper, const Params& params,
      bool throw_if_wrong_type = !std::is_same<Params, empty_pod>::value);
  any make_core_eval_pack(Operator oper) const;
 
  //-------
  // profiling
  //-------
  static const bool skip_core_ints = false;
};  // struct Engine
 
 
}  // namespace libint2
 
#ifndef LIBINT2_DOES_NOT_INLINE_ENGINE
#include "./engine.impl.h"
#endif
 
#endif /* _libint2_src_lib_libint_engine_h_ */