LIBINT  2.6.0
engine.h
1 /*
2  * Copyright (C) 2004-2019 Edward F. Valeev
3  *
4  * This file is part of Libint.
5  *
6  * Libint is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU Lesser General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * Libint is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public License
17  * along with Libint. If not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef _libint2_src_lib_libint_engine_h_
22 #define _libint2_src_lib_libint_engine_h_
23 
24 #ifndef LIBINT2_DOES_NOT_INLINE_ENGINE
25 # define __libint2_engine_inline inline
26 #else
27 # define __libint2_engine_inline
28 #endif
29 
30 #include <libint2/util/cxxstd.h>
31 #if LIBINT2_CPLUSPLUS_STD < 2011
32 # error "libint2/engine.h requires C++11 support"
33 #endif
34 
35 #include <algorithm>
36 #include <array>
37 #include <cstring>
38 #include <cstdio>
39 #include <functional>
40 #include <iostream>
41 #include <limits>
42 #include <map>
43 #include <memory>
44 #include <tuple>
45 #include <utility>
46 #include <vector>
47 
48 #include <libint2/cxxapi.h>
49 #include <libint2/boys_fwd.h>
50 #include <libint2/shell.h>
51 #include <libint2/solidharmonics.h>
52 #include <libint2/cartesian.h>
53 #include <libint2/util/any.h>
54 #include <libint2/util/array_adaptor.h>
55 #include <libint2/util/intpart_iter.h>
56 #include <libint2/util/compressed_pair.h>
57 #include <libint2/util/timer.h>
58 #include <libint2/cgshellinfo.h>
59 
60 // the engine will be profiled by default if library was configured with
61 // --enable-profile
62 #ifdef LIBINT2_PROFILE
63 #define LIBINT2_ENGINE_TIMERS
64 // uncomment if want to profile each integral class
65 #define LIBINT2_ENGINE_PROFILE_CLASS
66 #endif
67 // uncomment if want to profile the engine even if library was configured
68 // without --enable-profile
69 //# define LIBINT2_ENGINE_TIMERS
70 
71 namespace libint2 {
72 
76 typedef std::vector<std::pair<double, double>> ContractedGaussianGeminal;
77 
78 constexpr size_t num_geometrical_derivatives(size_t ncenter,
79  size_t deriv_order) {
80  return (deriv_order > 0)
81  ? (num_geometrical_derivatives(ncenter, deriv_order - 1) *
82  (3 * ncenter + deriv_order - 1)) /
83  deriv_order
84  : 1;
85 }
86 
87 template <typename T, unsigned N>
88 __libint2_engine_inline typename std::remove_all_extents<T>::type* to_ptr1(T (&a)[N]);
89 
91 
95 enum class Operator {
97  overlap = 0,
99  kinetic,
101  nuclear,
104  erf_nuclear,
107  erfc_nuclear,
112  emultipole1,
115  emultipole2,
118  emultipole3,
129  sphemultipole,
131  delta,
133  coulomb,
135  r12_m1 = coulomb,
137  cgtg,
139  cgtg_x_coulomb,
141  delcgtg2,
143  r12,
145  r12_1 = r12,
148  erf_coulomb,
151  erfc_coulomb,
153  stg,
155  stg_x_coulomb,
157  yukawa = stg_x_coulomb,
158  // do not modify this
159  invalid = -1,
160  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!keep this updated!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
161  first_1body_oper = overlap,
162  last_1body_oper = sphemultipole,
163  first_2body_oper = delta,
164  last_2body_oper = stg_x_coulomb,
165  first_oper = first_1body_oper,
166  last_oper = last_2body_oper
167 };
168 
172 inline constexpr int rank(Operator oper) {
173  return (oper >= Operator::first_1body_oper &&
174  oper <= Operator::last_1body_oper)
175  ? 1 :
176  ((oper >= Operator::first_2body_oper &&
177  oper <= Operator::last_2body_oper)
178  ? 2 : throw std::logic_error("rank(Operator): invalid operator given"));
179 }
180 
181 namespace detail {
182 struct default_operator_traits {
183  typedef struct {
184  } oper_params_type;
185  static oper_params_type default_params() { return oper_params_type{}; }
186  static constexpr auto nopers = 1u;
187  struct _core_eval_type {
188  template <typename... params>
189  static std::shared_ptr<const _core_eval_type> instance(params...) {
190  return nullptr;
191  }
192  };
193  using core_eval_type = const _core_eval_type;
194 };
195 } // namespace detail
196 
202 template <Operator Op>
203 struct operator_traits : public detail::default_operator_traits {};
204 
205 template <>
206 struct operator_traits<Operator::nuclear>
207  : public detail::default_operator_traits {
209  typedef std::vector<std::pair<scalar_type, std::array<scalar_type, 3>>>
210  oper_params_type;
211  static oper_params_type default_params() { return oper_params_type{}; }
212 #ifndef LIBINT_USER_DEFINED_REAL
213  typedef const libint2::FmEval_Chebyshev7<scalar_type> core_eval_type;
214 #else
215  typedef const libint2::FmEval_Reference<scalar_type> core_eval_type;
216 #endif
217 };
218 
219 template <>
220 struct operator_traits<Operator::erf_nuclear>
221  : public detail::default_operator_traits {
224  typedef std::tuple<
225  scalar_type, typename operator_traits<Operator::nuclear>::oper_params_type>
226  oper_params_type;
227  static oper_params_type default_params() {
228  return std::make_tuple(0,operator_traits<Operator::nuclear>::default_params());
229  }
231  core_eval_type;
232 };
233 
234 template <>
235 struct operator_traits<Operator::erfc_nuclear>
236  : public detail::default_operator_traits {
239  typedef typename operator_traits<Operator::erf_nuclear>::oper_params_type oper_params_type;
240  static oper_params_type default_params() {
241  return std::make_tuple(0,operator_traits<Operator::nuclear>::default_params());
242  }
244  core_eval_type;
245 };
246 
247 template <>
248 struct operator_traits<Operator::emultipole1>
249  : public detail::default_operator_traits {
252  typedef std::array<double, 3> oper_params_type;
253  static oper_params_type default_params() {
254  return oper_params_type{{0,0,0}};
255  }
256  static constexpr auto nopers = 4u;
257 };
258 template <>
259 struct operator_traits<Operator::emultipole2>
260  : public operator_traits<Operator::emultipole1> {
261  static constexpr auto nopers =
262  operator_traits<Operator::emultipole1>::nopers +
263  6;
264 };
265 template <>
266 struct operator_traits<Operator::emultipole3>
267  : public operator_traits<Operator::emultipole1> {
268  static constexpr auto nopers =
269  operator_traits<Operator::emultipole2>::nopers + 10;
270 };
271 template <>
272 struct operator_traits<Operator::sphemultipole>
273  : public operator_traits<Operator::emultipole1> {
274  static constexpr auto nopers =
275  (MULTIPOLE_MAX_ORDER + 1) * (MULTIPOLE_MAX_ORDER + 1);
276 };
277 
278 template <>
279 struct operator_traits<Operator::coulomb>
280  : public detail::default_operator_traits {
281 #ifndef LIBINT_USER_DEFINED_REAL
282  typedef const libint2::FmEval_Chebyshev7<scalar_type> core_eval_type;
283 #else
284  typedef const libint2::FmEval_Reference<scalar_type> core_eval_type;
285 #endif
286 };
287 namespace detail {
288 template <int K>
289 struct cgtg_operator_traits : public detail::default_operator_traits {
290  typedef libint2::GaussianGmEval<scalar_type, K> core_eval_type;
291  typedef ContractedGaussianGeminal oper_params_type;
292 };
293 } // namespace detail
294 template <>
295 struct operator_traits<Operator::cgtg>
296  : public detail::cgtg_operator_traits<0> {};
297 template <>
298 struct operator_traits<Operator::cgtg_x_coulomb>
299  : public detail::cgtg_operator_traits<-1> {};
300 template <>
301 struct operator_traits<Operator::delcgtg2>
302  : public detail::cgtg_operator_traits<2> {};
303 
304 template <>
305 struct operator_traits<Operator::delta>
306  : public detail::default_operator_traits {
308  core_eval_type;
309 };
310 
311 template <>
312 struct operator_traits<Operator::r12>
313  : public detail::default_operator_traits {
315  core_eval_type;
316 };
317 
318 template <>
319 struct operator_traits<Operator::erf_coulomb>
320  : public detail::default_operator_traits {
322  typedef scalar_type oper_params_type;
323  static oper_params_type default_params() {
324  return oper_params_type{0};
325  }
327  core_eval_type;
328 };
329 template <>
330 struct operator_traits<Operator::erfc_coulomb>
331  : public detail::default_operator_traits {
333  typedef scalar_type oper_params_type;
334  static oper_params_type default_params() {
335  return oper_params_type{0};
336  }
338  core_eval_type;
339 };
340 
341 template <>
342 struct operator_traits<Operator::stg>
343  : public detail::default_operator_traits {
345  typedef scalar_type oper_params_type;
346  static oper_params_type default_params() {
347  return oper_params_type{0};
348  }
350  core_eval_type;
351 };
352 
353 template <>
354 struct operator_traits<Operator::stg_x_coulomb>
355  : public detail::default_operator_traits {
357  typedef scalar_type oper_params_type;
358  static oper_params_type default_params() {
359  return oper_params_type{0};
360  }
362  core_eval_type;
363 };
364 
367 default_params(const Operator& oper);
368 
369 namespace detail {
370 template <typename core_eval_type>
371 using __core_eval_pack_type =
372  compressed_pair<std::shared_ptr<core_eval_type>,
374 template <Operator Op>
375 using core_eval_pack_type =
376  __core_eval_pack_type<typename operator_traits<Op>::core_eval_type>;
377 }
378 
385 enum class BraKet {
386  x_x = 0,
387  xx_xx,
388  xs_xx,
389  xx_xs,
390  xs_xs,
391  invalid = -1,
392  first_1body_braket = x_x,
393  last_1body_braket = x_x,
394  first_2body_braket = xx_xx,
395  last_2body_braket = xs_xs,
396  first_braket = first_1body_braket,
397  last_braket = last_2body_braket
398 };
399 #define BOOST_PP_NBODY_BRAKET_MAX_INDEX 4
400 
404 inline constexpr int rank(BraKet braket) {
405  return (braket==BraKet::x_x || braket==BraKet::xs_xs) ? 2 :
406  ((braket==BraKet::xs_xx || braket==BraKet::xx_xs) ? 3 :
407  ((braket==BraKet::xx_xx) ? 4 : throw std::logic_error("rank(BraKet): invalid braket given")));
408 }
409 
413 inline constexpr BraKet default_braket(Operator oper) {
414  return (rank(oper)==1) ? BraKet::x_x :
415  ((rank(oper)==2) ? BraKet::xx_xx : throw std::logic_error("default_braket(Operator): invalid operator given"));
416 }
417 
418 constexpr size_t nopers_2body = static_cast<int>(Operator::last_2body_oper) -
419  static_cast<int>(Operator::first_2body_oper) +
420  1;
421 constexpr size_t nbrakets_2body = static_cast<int>(BraKet::last_2body_braket) -
422  static_cast<int>(BraKet::first_2body_braket) +
423  1;
424 constexpr size_t nderivorders_2body = LIBINT2_MAX_DERIV_ORDER + 1;
425 
431 class Engine {
432  private:
433  typedef struct {
434  } empty_pod;
435 
436  public:
437  static constexpr auto max_ntargets =
438  std::extent<decltype(std::declval<Libint_t>().targets), 0>::value;
439  using target_ptr_vec =
440  std::vector<const value_type*, detail::ext_stack_allocator<const value_type*, max_ntargets>>;
441 
446  Engine()
447  : oper_(Operator::invalid),
448  braket_(BraKet::invalid),
449  primdata_(),
450  stack_size_(0),
451  lmax_(-1),
452  deriv_order_(0),
453  cartesian_shell_normalization_(CartesianShellNormalization::standard) {
454  set_precision(std::numeric_limits<scalar_type>::epsilon());
455  }
456 
457  // clang-format off
459 
489  // clang-format on
490  template <typename Params = empty_pod>
491  Engine(Operator oper, size_t max_nprim, int max_l, int deriv_order = 0,
492  scalar_type precision = std::numeric_limits<scalar_type>::epsilon(),
493  Params params = empty_pod(), BraKet braket = BraKet::invalid)
494  : oper_(oper),
495  braket_(braket),
496  primdata_(),
497  spbra_(max_nprim),
498  spket_(max_nprim),
499  stack_size_(0),
500  lmax_(max_l),
501  deriv_order_(deriv_order),
502  cartesian_shell_normalization_(CartesianShellNormalization::standard),
503  params_(enforce_params_type(oper, params)) {
504  set_precision(precision);
505  assert(max_nprim > 0);
506  initialize(max_nprim);
507  core_eval_pack_ = make_core_eval_pack(oper); // must follow initialize() to
508  // ensure default braket_ has
509  // been set
510  init_core_ints_params(params_);
511  }
512 
514  // intel does not accept "move ctor = default"
515  Engine(Engine&& other)
516  : oper_(other.oper_),
517  braket_(other.braket_),
518  primdata_(std::move(other.primdata_)),
519  spbra_(std::move(other.spbra_)),
520  spket_(std::move(other.spket_)),
521  stack_size_(other.stack_size_),
522  lmax_(other.lmax_),
523  hard_lmax_(other.hard_lmax_),
524  hard_default_lmax_(other.hard_default_lmax_),
525  deriv_order_(other.deriv_order_),
526  precision_(other.precision_),
527  ln_precision_(other.ln_precision_),
528  cartesian_shell_normalization_(other.cartesian_shell_normalization_),
529  core_eval_pack_(std::move(other.core_eval_pack_)),
530  params_(std::move(other.params_)),
531  core_ints_params_(std::move(other.core_ints_params_)),
532  targets_(std::move(other.targets_)),
533  set_targets_(other.set_targets_),
534  scratch_(std::move(other.scratch_)),
535  scratch2_(other.scratch2_),
536  buildfnptrs_(other.buildfnptrs_) {}
537 
539  Engine(const Engine& other)
540  : oper_(other.oper_),
541  braket_(other.braket_),
542  primdata_(other.primdata_.size()),
543  spbra_(other.spbra_),
544  spket_(other.spket_),
545  stack_size_(other.stack_size_),
546  lmax_(other.lmax_),
547  deriv_order_(other.deriv_order_),
548  precision_(other.precision_),
549  ln_precision_(other.ln_precision_),
550  cartesian_shell_normalization_(other.cartesian_shell_normalization_),
551  core_eval_pack_(other.core_eval_pack_),
552  params_(other.params_),
553  core_ints_params_(other.core_ints_params_) {
554  initialize();
555  }
556 
557  ~Engine() { finalize(); }
558 
560  Engine& operator=(Engine&& other) {
561  oper_ = other.oper_;
562  braket_ = other.braket_;
563  primdata_ = std::move(other.primdata_);
564  spbra_ = std::move(other.spbra_);
565  spket_ = std::move(other.spket_);
566  stack_size_ = other.stack_size_;
567  lmax_ = other.lmax_;
568  hard_lmax_ = other.hard_lmax_;
569  hard_default_lmax_ = other.hard_default_lmax_;
570  deriv_order_ = other.deriv_order_;
571  precision_ = other.precision_;
572  ln_precision_ = other.ln_precision_;
573  cartesian_shell_normalization_ = other.cartesian_shell_normalization_;
574  core_eval_pack_ = std::move(other.core_eval_pack_);
575  params_ = std::move(other.params_);
576  core_ints_params_ = std::move(other.core_ints_params_);
577  targets_ = std::move(other.targets_);
578  set_targets_ = other.set_targets_;
579  scratch_ = std::move(other.scratch_);
580  scratch2_ = other.scratch2_;
581  buildfnptrs_ = other.buildfnptrs_;
582  return *this;
583  }
584 
586  Engine& operator=(const Engine& other) {
587  oper_ = other.oper_;
588  braket_ = other.braket_;
589  primdata_.resize(other.primdata_.size());
590  spbra_ = other.spbra_;
591  spket_ = other.spket_;
592  stack_size_ = other.stack_size_;
593  lmax_ = other.lmax_;
594  deriv_order_ = other.deriv_order_;
595  precision_ = other.precision_;
596  ln_precision_ = other.ln_precision_;
597  cartesian_shell_normalization_ = other.cartesian_shell_normalization_;
598  core_eval_pack_ = other.core_eval_pack_;
599  params_ = other.params_;
600  core_ints_params_ = other.core_ints_params_;
601  initialize();
602  return *this;
603  }
604 
606  int operator_rank() const { return rank(oper_); }
607 
609  int braket_rank() const { return rank(braket_); }
610 
612  Operator oper() const { return oper_; }
613 
615  BraKet braket() const { return braket_; }
616 
620  Engine& set(Operator new_oper) {
621  if (oper_ != new_oper) {
622  oper_ = new_oper;
623  braket_ = default_braket(oper_);
624  params_ = default_params(oper_);
625  initialize();
626  core_eval_pack_ = make_core_eval_pack(oper_); // must follow initialize() to
627  // ensure default braket_ has
628  // been set
629  }
630  return *this;
631  }
632 
635  Engine& set(BraKet new_braket) {
636  if (braket_ != new_braket) {
637  braket_ = new_braket;
638  initialize();
639  }
640  return *this;
641  }
642 
646  template <typename Params>
647  Engine& set_params(const Params& params) {
648  params_ = params;
649  init_core_ints_params(params_);
650  reset_scratch();
651  return *this;
652  }
653 
655  std::size_t max_nprim() const {
656  assert(spbra_.primpairs.size() == spket_.primpairs.size());
657  return static_cast<std::size_t>(std::sqrt(spbra_.primpairs.size()));
658  }
659 
663  Engine& set_max_nprim(std::size_t n) {
664  if (n*n > spbra_.primpairs.size()) {
665  spbra_.resize(n);
666  spket_.resize(n);
667  initialize(n);
668  }
669  return *this;
670  }
671 
673  std::size_t max_l() const {
674  return lmax_;
675  }
676 
680  Engine& set_max_l(std::size_t L) {
681  if (L >= static_cast<std::size_t>(lmax_)) {
682  lmax_ = L;
683  initialize();
684  }
685  return *this;
686  }
687 
693  const target_ptr_vec& results() const { return targets_; }
694 
698  unsigned int nshellsets() const { return targets_.size(); }
699 
701 
706  template <typename... ShellPack>
707  __libint2_engine_inline const target_ptr_vec& compute(
708  const libint2::Shell& first_shell, const ShellPack&... rest_of_shells);
709 
715  __libint2_engine_inline const target_ptr_vec& compute1(
716  const libint2::Shell& s1, const libint2::Shell& s2);
717 
739  template <Operator oper, BraKet braket, size_t deriv_order>
740  __libint2_engine_inline const target_ptr_vec& compute2(const Shell& bra1,
741  const Shell& bra2,
742  const Shell& ket1,
743  const Shell& ket2,
744  const ShellPair* spbra = nullptr,
745  const ShellPair* spket = nullptr);
746 
747  typedef const target_ptr_vec& (Engine::*compute2_ptr_type)(const Shell& bra1,
748  const Shell& bra2,
749  const Shell& ket1,
750  const Shell& ket2,
751  const ShellPair* spbra,
752  const ShellPair* spket);
753 
766  Engine& set_precision(scalar_type prec) {
767  if (prec <= 0.) {
768  precision_ = 0.;
769  ln_precision_ = std::numeric_limits<scalar_type>::lowest();
770  } else {
771  precision_ = prec;
772  using std::log;
773  ln_precision_ = log(precision_);
774  }
775  return *this;
776  }
779  scalar_type precision() const { return precision_; }
780 
782  CartesianShellNormalization cartesian_shell_normalization() const {
783  return cartesian_shell_normalization_;
784  }
785 
789  Engine& set(CartesianShellNormalization norm) {
790  cartesian_shell_normalization_ = norm;
791  return *this;
792  }
793 
795  void print_timers(std::ostream& os = std::cout) {
796 #ifdef LIBINT2_ENGINE_TIMERS
797  os << "timers: prereq = " << timers.read(0);
798 #ifndef LIBINT2_PROFILE // if libint's profiling was on, engine's build timer
799  // will include its overhead
800  // do not report it, detailed profiling from libint will be printed below
801  os << " build = " << timers.read(1);
802 #endif
803  os << " tform = " << timers.read(2) << std::endl;
804 #endif
805 #ifdef LIBINT2_PROFILE
806  os << "build timers: hrr = " << primdata_[0].timers->read(0)
807  << " vrr = " << primdata_[0].timers->read(1) << std::endl;
808 #endif
809 #ifdef LIBINT2_ENGINE_TIMERS
810 #ifdef LIBINT2_ENGINE_PROFILE_CLASS
811  for (const auto& p : class_profiles) {
812  char buf[1024];
813  std::snprintf(buf, sizeof(buf), "{\"%s\", %10.5lf, %10.5lf, %10.5lf, %10.5lf, %ld, %ld},",
814  p.first.to_string().c_str(), p.second.prereqs, p.second.build_vrr,
815  p.second.build_hrr, p.second.tform, p.second.nshellset,
816  p.second.nprimset);
817  os << buf << std::endl;
818  }
819 #endif
820 #endif
821  }
822 
824  class lmax_exceeded : virtual public std::logic_error {
825  public:
826  lmax_exceeded(const char* task_name, size_t lmax_limit,
827  size_t lmax_requested)
828  : std::logic_error(
829  "Engine::lmax_exceeded -- angular momentum limit exceeded"),
830  lmax_limit_(lmax_limit),
831  lmax_requested_(lmax_requested) {
832  strncpy(task_name_, task_name, 64);
833  task_name_[64] = '\0';
834  }
835  ~lmax_exceeded() noexcept {}
836 
837  const char* task_name() const {
838  return static_cast<const char*>(task_name_);
839  }
840  size_t lmax_limit() const { return lmax_limit_; }
841  size_t lmax_requested() const { return lmax_requested_; }
842 
843  private:
844  char task_name_[65];
845  size_t lmax_limit_;
846  size_t lmax_requested_;
847  };
848 
850  class using_default_initialized : virtual public std::logic_error {
851  public:
852  using_default_initialized()
853  : std::logic_error(
854  "Engine::using_default_initialized -- attempt to use a default-initialized Engine") {
855  }
856  ~using_default_initialized() noexcept {}
857  };
858 
859  private:
860  Operator oper_;
861  BraKet braket_;
862  std::vector<Libint_t> primdata_;
863  ShellPair spbra_, spket_;
864  size_t stack_size_; // amount allocated by libint2_init_xxx in
865  // primdata_[0].stack
866  int lmax_;
867  int hard_lmax_; // max L supported by library for this operator type (i.e. LIBINT2_MAX_AM_<task><deriv_order>) + 1
868  int hard_default_lmax_; // max L supported by library for default operator type
869  // (i.e. LIBINT2_MAX_AM_default<deriv_order>) + 1
870  // this is only set and used if LIBINT2_CENTER_DEPENDENT_MAX_AM_<task><deriv_order> == 1
871  int deriv_order_;
872  scalar_type precision_;
873  scalar_type ln_precision_;
874 
875  // specifies the normalization convention for Cartesian Gaussians
876  CartesianShellNormalization cartesian_shell_normalization_;
877 
878  any core_eval_pack_;
879 
881  any params_; // operator params
886  any core_ints_params_;
888  void init_core_ints_params(const any& params);
889 
892  target_ptr_vec targets_;
895  bool set_targets_;
896 
897  std::vector<value_type>
898  scratch_; // for transposes and/or transforming to solid harmonics
899  value_type* scratch2_; // &scratch_[0] points to the first block large enough to
900  // hold all target ints
901  // scratch2_ points to second such block. It could point into scratch_ or at
902  // primdata_[0].stack
903  typedef void (*buildfnptr_t)(const Libint_t*);
904  buildfnptr_t* buildfnptrs_;
905 
907  unsigned int compute_nshellsets() const {
908  const unsigned int num_operator_geometrical_derivatives =
909  (oper_ == Operator::nuclear || oper_ == Operator::erf_nuclear ||
910  oper_ == Operator::erfc_nuclear)
911  ? this->nparams()
912  : 0;
913  const auto ncenters = braket_rank() + num_operator_geometrical_derivatives;
914  return nopers() * num_geometrical_derivatives(ncenters, deriv_order_);
915  }
916 
917  void reset_scratch() {
918  const auto nshsets = compute_nshellsets();
919  targets_.resize(nshsets);
920  set_targets_ = (&targets_[0] != const_cast<const value_type**>(primdata_[0].targets));
921  const auto ncart_max = (lmax_ + 1) * (lmax_ + 2) / 2;
922  const auto target_shellset_size =
923  nshsets * std::pow(ncart_max, braket_rank());
924  // need to be able to hold 2 sets of target shellsets: the worst case occurs
925  // when dealing with
926  // 1-body Coulomb ints derivatives ... have 2+natom 1st-order derivative sets
927  // (and many more of 2nd and higher) that are stored in scratch
928  // then need to transform to solids. To avoid copying back and forth make
929  // sure that there is enough
930  // room to transform all ints and save them in correct order in single pass
931  const auto need_extra_large_scratch = stack_size_ < target_shellset_size;
932  scratch_.resize(need_extra_large_scratch ? 2 * target_shellset_size
933  : target_shellset_size);
934  scratch2_ = need_extra_large_scratch ? &scratch_[target_shellset_size]
935  : primdata_[0].stack;
936  }
937 
938  __libint2_engine_inline void compute_primdata(Libint_t& primdata,
939  const Shell& s1,
940  const Shell& s2, size_t p1,
941  size_t p2, size_t oset);
942 
945  __libint2_engine_inline const std::vector<Engine::compute2_ptr_type>&
946  compute2_ptrs() const;
947 
948  // max_nprim=0 avoids resizing primdata_
949  __libint2_engine_inline void initialize(size_t max_nprim = 0);
950  // generic _initializer
951  __libint2_engine_inline void _initialize();
952 
953  void finalize() {
954  if (primdata_.size() != 0) {
955  libint2_cleanup_default(&primdata_[0]);
956  }
957  } // finalize()
958 
959  //-------
960  // utils
961  //-------
962  unsigned int nparams() const;
963  unsigned int nopers() const;
970  template <typename Params>
971  static any enforce_params_type(
972  Operator oper, const Params& params,
973  bool throw_if_wrong_type = !std::is_same<Params, empty_pod>::value);
976  any make_core_eval_pack(Operator oper) const;
977 
978  //-------
979  // profiling
980  //-------
981  static const bool skip_core_ints = false;
982 }; // struct Engine
983 
984 
985 } // namespace libint2
986 
987 #ifndef LIBINT2_DOES_NOT_INLINE_ENGINE
988 #include "./engine.impl.h"
989 #endif
990 
991 #endif /* _libint2_src_lib_libint_engine_h_ */
992 
Definition: boys_fwd.h:50
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition: shell.h:81
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
some evaluators need thread-local scratch, but most don't
Definition: boys.h:1331
CartesianShellNormalization
Normalization convention for Cartesian Gaussian shells.
Definition: cgshellinfo.h:31
Computes the Boys function, $ F_m (T) = \int_0^1 u^{2m} \exp(-T u^2) \, {\rm d}u $,...
Definition: boys.h:247
void initialize(bool verbose=false)
initializes the libint library if not currently initialized, no-op otherwise
Definition: initialize.h:68
Definition: boys.h:1326
Array idential to C++0X arrays.
Definition: stdarray_bits.h:14
Computes the Boys function, , using single algorithm (asymptotic expansion).
Definition: boys.h:144
core integral for Yukawa and exponential interactions
Definition: boys.h:834
void finalize()
finalizes the libint library.
Definition: initialize.h:88
a partial C++17 std::any implementation (and less efficient than can be)
Definition: any.h:67
__libint2_engine_inline libint2::any default_params(const Operator &oper)
the runtime version of operator_traits<oper>::default_params()
Definition: engine.impl.h:110