21 #ifndef _libint2_src_lib_libint_engine_h_ 22 #define _libint2_src_lib_libint_engine_h_ 24 #ifndef LIBINT2_DOES_NOT_INLINE_ENGINE 25 # define __libint2_engine_inline inline 27 # define __libint2_engine_inline 30 #include <libint2/util/cxxstd.h> 31 #if LIBINT2_CPLUSPLUS_STD < 2011 32 # error "libint2/engine.h requires C++11 support" 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> 62 #ifdef LIBINT2_PROFILE 63 #define LIBINT2_ENGINE_TIMERS 65 #define LIBINT2_ENGINE_PROFILE_CLASS 76 typedef std::vector<std::pair<double, double>> ContractedGaussianGeminal;
78 constexpr
size_t num_geometrical_derivatives(
size_t ncenter,
80 return (deriv_order > 0)
81 ? (num_geometrical_derivatives(ncenter, deriv_order - 1) *
82 (3 * ncenter + deriv_order - 1)) /
87 template <
typename T,
unsigned N>
88 __libint2_engine_inline
typename std::remove_all_extents<T>::type* to_ptr1(T (&a)[N]);
157 yukawa = stg_x_coulomb,
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
172 inline constexpr
int rank(Operator oper) {
173 return (oper >= Operator::first_1body_oper &&
174 oper <= Operator::last_1body_oper)
176 ((oper >= Operator::first_2body_oper &&
177 oper <= Operator::last_2body_oper)
178 ? 2 :
throw std::logic_error(
"rank(Operator): invalid operator given"));
182 struct default_operator_traits {
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...) {
193 using core_eval_type =
const _core_eval_type;
202 template <Operator Op>
203 struct operator_traits :
public detail::default_operator_traits {};
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>>>
211 static oper_params_type
default_params() {
return oper_params_type{}; }
212 #ifndef LIBINT_USER_DEFINED_REAL 220 struct operator_traits<Operator::erf_nuclear>
221 :
public detail::default_operator_traits {
225 scalar_type,
typename operator_traits<Operator::nuclear>::oper_params_type>
228 return std::make_tuple(0,operator_traits<Operator::nuclear>::default_params());
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;
241 return std::make_tuple(0,operator_traits<Operator::nuclear>::default_params());
248 struct operator_traits<Operator::emultipole1>
249 :
public detail::default_operator_traits {
254 return oper_params_type{{0,0,0}};
256 static constexpr
auto nopers = 4u;
259 struct operator_traits<Operator::emultipole2>
260 :
public operator_traits<Operator::emultipole1> {
261 static constexpr
auto nopers =
262 operator_traits<Operator::emultipole1>::nopers +
266 struct operator_traits<Operator::emultipole3>
267 :
public operator_traits<Operator::emultipole1> {
268 static constexpr
auto nopers =
269 operator_traits<Operator::emultipole2>::nopers + 10;
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);
279 struct operator_traits<Operator::coulomb>
280 :
public detail::default_operator_traits {
281 #ifndef LIBINT_USER_DEFINED_REAL 289 struct cgtg_operator_traits :
public detail::default_operator_traits {
291 typedef ContractedGaussianGeminal oper_params_type;
295 struct operator_traits<Operator::cgtg>
296 :
public detail::cgtg_operator_traits<0> {};
298 struct operator_traits<Operator::cgtg_x_coulomb>
299 :
public detail::cgtg_operator_traits<-1> {};
301 struct operator_traits<Operator::delcgtg2>
302 :
public detail::cgtg_operator_traits<2> {};
305 struct operator_traits<Operator::delta>
306 :
public detail::default_operator_traits {
312 struct operator_traits<Operator::r12>
313 :
public detail::default_operator_traits {
319 struct operator_traits<Operator::erf_coulomb>
320 :
public detail::default_operator_traits {
322 typedef scalar_type oper_params_type;
324 return oper_params_type{0};
330 struct operator_traits<Operator::erfc_coulomb>
331 :
public detail::default_operator_traits {
333 typedef scalar_type oper_params_type;
335 return oper_params_type{0};
342 struct operator_traits<Operator::stg>
343 :
public detail::default_operator_traits {
345 typedef scalar_type oper_params_type;
347 return oper_params_type{0};
354 struct operator_traits<Operator::stg_x_coulomb>
355 :
public detail::default_operator_traits {
357 typedef scalar_type oper_params_type;
359 return oper_params_type{0};
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>;
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
399 #define BOOST_PP_NBODY_BRAKET_MAX_INDEX 4 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")));
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"));
418 constexpr
size_t nopers_2body = static_cast<int>(Operator::last_2body_oper) -
419 static_cast<int>(Operator::first_2body_oper) +
421 constexpr
size_t nbrakets_2body = static_cast<int>(BraKet::last_2body_braket) -
422 static_cast<int>(BraKet::first_2body_braket) +
424 constexpr
size_t nderivorders_2body = LIBINT2_MAX_DERIV_ORDER + 1;
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>>;
447 : oper_(Operator::invalid),
448 braket_(BraKet::invalid),
454 set_precision(std::numeric_limits<scalar_type>::epsilon());
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)
501 deriv_order_(deriv_order),
503 params_(enforce_params_type(oper, params)) {
504 set_precision(precision);
505 assert(max_nprim > 0);
507 core_eval_pack_ = make_core_eval_pack(oper);
510 init_core_ints_params(params_);
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_),
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_) {}
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_),
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_) {
560 Engine& operator=(Engine&& other) {
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_;
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_;
586 Engine& operator=(
const Engine& other) {
588 braket_ = other.braket_;
589 primdata_.resize(other.primdata_.size());
590 spbra_ = other.spbra_;
591 spket_ = other.spket_;
592 stack_size_ = other.stack_size_;
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_;
606 int operator_rank()
const {
return rank(oper_); }
609 int braket_rank()
const {
return rank(braket_); }
612 Operator oper()
const {
return oper_; }
615 BraKet braket()
const {
return braket_; }
620 Engine& set(Operator new_oper) {
621 if (oper_ != new_oper) {
623 braket_ = default_braket(oper_);
626 core_eval_pack_ = make_core_eval_pack(oper_);
635 Engine& set(BraKet new_braket) {
636 if (braket_ != new_braket) {
637 braket_ = new_braket;
646 template <
typename Params>
647 Engine& set_params(
const Params& params) {
649 init_core_ints_params(params_);
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()));
663 Engine& set_max_nprim(std::size_t n) {
664 if (n*n > spbra_.primpairs.size()) {
673 std::size_t max_l()
const {
680 Engine& set_max_l(std::size_t L) {
681 if (L >= static_cast<std::size_t>(lmax_)) {
693 const target_ptr_vec& results()
const {
return targets_; }
698 unsigned int nshellsets()
const {
return targets_.size(); }
706 template <
typename... ShellPack>
707 __libint2_engine_inline
const target_ptr_vec& compute(
708 const libint2::Shell& first_shell,
const ShellPack&... rest_of_shells);
715 __libint2_engine_inline
const target_ptr_vec& compute1(
739 template <Operator oper, BraKet braket,
size_t deriv_order>
740 __libint2_engine_inline
const target_ptr_vec& compute2(
const Shell& bra1,
744 const ShellPair* spbra =
nullptr,
745 const ShellPair* spket =
nullptr);
747 typedef const target_ptr_vec& (Engine::*compute2_ptr_type)(
const Shell& bra1,
751 const ShellPair* spbra,
752 const ShellPair* spket);
766 Engine& set_precision(scalar_type prec) {
769 ln_precision_ = std::numeric_limits<scalar_type>::lowest();
773 ln_precision_ = log(precision_);
779 scalar_type precision()
const {
return precision_; }
783 return cartesian_shell_normalization_;
790 cartesian_shell_normalization_ = norm;
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 801 os <<
" build = " << timers.read(1);
803 os <<
" tform = " << timers.read(2) << std::endl;
805 #ifdef LIBINT2_PROFILE 806 os <<
"build timers: hrr = " << primdata_[0].timers->read(0)
807 <<
" vrr = " << primdata_[0].timers->read(1) << std::endl;
809 #ifdef LIBINT2_ENGINE_TIMERS 810 #ifdef LIBINT2_ENGINE_PROFILE_CLASS 811 for (
const auto& p : class_profiles) {
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,
817 os << buf << std::endl;
824 class lmax_exceeded :
virtual public std::logic_error {
826 lmax_exceeded(
const char* task_name,
size_t lmax_limit,
827 size_t lmax_requested)
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';
835 ~lmax_exceeded() noexcept {}
837 const char* task_name()
const {
838 return static_cast<const char*>(task_name_);
840 size_t lmax_limit()
const {
return lmax_limit_; }
841 size_t lmax_requested()
const {
return lmax_requested_; }
846 size_t lmax_requested_;
850 class using_default_initialized :
virtual public std::logic_error {
852 using_default_initialized()
854 "Engine::using_default_initialized -- attempt to use a default-initialized Engine") {
856 ~using_default_initialized() noexcept {}
862 std::vector<Libint_t> primdata_;
863 ShellPair spbra_, spket_;
868 int hard_default_lmax_;
872 scalar_type precision_;
873 scalar_type ln_precision_;
886 any core_ints_params_;
888 void init_core_ints_params(
const any& params);
892 target_ptr_vec targets_;
897 std::vector<value_type>
899 value_type* scratch2_;
903 typedef void (*buildfnptr_t)(
const Libint_t*);
904 buildfnptr_t* buildfnptrs_;
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)
913 const auto ncenters = braket_rank() + num_operator_geometrical_derivatives;
914 return nopers() * num_geometrical_derivatives(ncenters, deriv_order_);
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());
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;
938 __libint2_engine_inline
void compute_primdata(Libint_t& primdata,
940 const Shell& s2,
size_t p1,
941 size_t p2,
size_t oset);
945 __libint2_engine_inline
const std::vector<Engine::compute2_ptr_type>&
946 compute2_ptrs()
const;
949 __libint2_engine_inline
void initialize(
size_t max_nprim = 0);
951 __libint2_engine_inline
void _initialize();
954 if (primdata_.size() != 0) {
955 libint2_cleanup_default(&primdata_[0]);
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;
981 static const bool skip_core_ints =
false;
987 #ifndef LIBINT2_DOES_NOT_INLINE_ENGINE 988 #include "./engine.impl.h" 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
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