21 #ifndef _libint2_src_lib_libint_engineimpl_h_ 22 #define _libint2_src_lib_libint_engineimpl_h_ 28 #pragma GCC diagnostic push 29 #pragma GCC system_header 31 #pragma GCC diagnostic pop 33 #include <libint2/boys.h> 34 #if LIBINT_HAS_SYSTEM_BOOST_PREPROCESSOR_VARIADICS 35 # include <boost/preprocessor.hpp> 36 # include <boost/preprocessor/facilities/is_1.hpp> 37 #else // use bundled boost 38 # include <libint2/boost/preprocessor.hpp> 39 # include <libint2/boost/preprocessor/facilities/is_1.hpp> 44 #define BOOST_PP_MAKE_TUPLE_INTERNAL(z, i, last) \ 45 i BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, last)) 46 #define BOOST_PP_MAKE_TUPLE(n) \ 48 (BOOST_PP_REPEAT(n, BOOST_PP_MAKE_TUPLE_INTERNAL, BOOST_PP_DEC(n))) 52 #ifdef LIBINT2_PROFILE 53 #define LIBINT2_ENGINE_TIMERS 55 #define LIBINT2_ENGINE_PROFILE_CLASS 63 template <
typename T,
unsigned N>
64 typename std::remove_all_extents<T>::type* to_ptr1(T (&a)[N]) {
65 return reinterpret_cast<typename std::remove_all_extents<T>::type*
>(&a);
72 #define BOOST_PP_NBODY_OPERATOR_LIST \ 82 (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, (eri, BOOST_PP_NIL))))))))))))))))))) 84 #define BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE \ 85 BOOST_PP_MAKE_TUPLE(BOOST_PP_LIST_SIZE(BOOST_PP_NBODY_OPERATOR_LIST)) 86 #define BOOST_PP_NBODY_OPERATOR_INDEX_LIST \ 87 BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE) 88 #define BOOST_PP_NBODY_OPERATOR_LAST_ONEBODY_INDEX \ 89 8 // sphemultipole, the 9th member of BOOST_PP_NBODY_OPERATOR_LIST, is the last 93 #define BOOST_PP_NBODY_BRAKET_INDEX_TUPLE \ 94 BOOST_PP_MAKE_TUPLE(BOOST_PP_INC(BOOST_PP_NBODY_BRAKET_MAX_INDEX)) 95 #define BOOST_PP_NBODY_BRAKET_INDEX_LIST \ 96 BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_BRAKET_INDEX_TUPLE) 97 #define BOOST_PP_NBODY_BRAKET_RANK_TUPLE (2, 3, 4) 98 #define BOOST_PP_NBODY_BRAKET_RANK_LIST \ 99 BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_BRAKET_RANK_TUPLE) 102 #define BOOST_PP_NBODY_DERIV_ORDER_TUPLE \ 103 BOOST_PP_MAKE_TUPLE(BOOST_PP_INC(LIBINT2_MAX_DERIV_ORDER)) 104 #define BOOST_PP_NBODY_DERIV_ORDER_LIST \ 105 BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_DERIV_ORDER_TUPLE) 111 switch (static_cast<int>(oper)) {
112 #define BOOST_PP_NBODYENGINE_MCR1(r, data, i, elem) \ 114 return operator_traits<static_cast<Operator>(i)>::default_params(); 115 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR1, _,
116 BOOST_PP_NBODY_OPERATOR_LIST)
120 assert(
false &&
"missing case in switch");
133 template <
typename... ShellPack>
134 __libint2_engine_inline
const Engine::target_ptr_vec& Engine::compute(
135 const libint2::Shell& first_shell,
const ShellPack&... rest_of_shells) {
136 constexpr
auto nargs = 1 +
sizeof...(rest_of_shells);
137 assert(nargs == braket_rank() &&
"# of arguments to compute() does not match the braket type");
140 first_shell, rest_of_shells...}};
142 if (operator_rank() == 1) {
143 if (nargs == 2)
return compute1(shells[0], shells[1]);
144 }
else if (operator_rank() == 2) {
145 auto compute_ptr_idx = ((static_cast<int>(oper_) -
146 static_cast<int>(Operator::first_2body_oper)) *
148 (static_cast<int>(braket_) -
149 static_cast<int>(BraKet::first_2body_braket))) *
152 auto compute_ptr = compute2_ptrs().at(compute_ptr_idx);
153 assert(compute_ptr !=
nullptr &&
"2-body compute function not found");
155 return (this->*compute_ptr)(shells[0],
Shell::unit(), shells[1],
158 return (this->*compute_ptr)(shells[0],
Shell::unit(), shells[1],
159 shells[2],
nullptr,
nullptr);
161 return (this->*compute_ptr)(shells[0], shells[1], shells[2], shells[3],
nullptr,
nullptr);
164 assert(
false &&
"missing feature");
172 __libint2_engine_inline
const Engine::target_ptr_vec& Engine::compute1(
175 assert((s1.ncontr() == 1 && s2.ncontr() == 1) &&
176 "generally-contracted shells not yet supported");
178 const auto oper_is_nuclear =
179 (oper_ == Operator::nuclear || oper_ == Operator::erf_nuclear ||
180 oper_ == Operator::erfc_nuclear);
182 const auto l1 = s1.
contr[0].l;
183 const auto l2 = s2.
contr[0].l;
184 assert(l1 <= lmax_ &&
"the angular momentum limit is exceeded");
185 assert(l2 <= lmax_ &&
"the angular momentum limit is exceeded");
189 if (oper_is_nuclear && nparams() == 0)
190 throw std::logic_error(
191 "Engine<*nuclear>, but no charges found; forgot to call " 194 const auto n1 = s1.size();
195 const auto n2 = s2.size();
196 const auto n12 = n1 * n2;
197 const auto ncart1 = s1.cartesian_size();
198 const auto ncart2 = s2.cartesian_size();
199 const auto ncart12 = ncart1 * ncart2;
202 const auto nprim1 = s1.nprim();
203 const auto nprim2 = s2.nprim();
204 const auto nprimpairs = nprim1 * nprim2;
205 assert(nprimpairs <= primdata_.size() &&
"the max number of primitive pairs exceeded");
207 auto nparam_sets = nparams();
210 bool set_targets = set_targets_;
213 const auto ntargets = nopers() * num_geometrical_derivatives(2, deriv_order_);
220 const auto nderivcenters_shset = 2 + (oper_is_nuclear ? nparam_sets : 0);
221 const auto nderivcoord = 3 * nderivcenters_shset;
222 const auto num_shellsets_computed =
223 nopers() * num_geometrical_derivatives(nderivcenters_shset, deriv_order_);
231 const auto accumulate_ints_in_scratch = oper_is_nuclear;
234 const auto lmax = std::max(l1, l2);
235 assert(lmax <= lmax_ &&
"the angular momentum limit is exceeded");
241 const auto tform_to_solids =
242 (s1.
contr[0].pure || s2.
contr[0].pure) && lmax != 0;
246 const auto compute_directly =
247 lmax == 0 && deriv_order_ == 0 &&
248 (oper_ == Operator::overlap || oper_is_nuclear);
249 if (compute_directly) {
250 primdata_[0].stack[0] = 0;
251 targets_[0] = primdata_[0].stack;
254 if (accumulate_ints_in_scratch)
255 std::fill(std::begin(scratch_),
256 std::begin(scratch_) + num_shellsets_computed * ncart12, 0.0);
259 for (
auto pset = 0u; pset != nparam_sets; ++pset) {
260 if (!oper_is_nuclear)
261 assert(nparam_sets == 1 &&
"unexpected number of operator parameters");
264 for (
auto p1 = 0; p1 != nprim1; ++p1) {
265 for (
auto p2 = 0; p2 != nprim2; ++p2, ++p12) {
266 compute_primdata(primdata_[p12], s1, s2, p1, p2, pset);
269 primdata_[0].contrdepth = p12;
271 if (compute_directly) {
272 auto& result = primdata_[0].stack[0];
274 case Operator::overlap:
275 for (
auto p12 = 0; p12 != primdata_[0].contrdepth; ++p12)
276 result += primdata_[p12]._0_Overlap_0_x[0] *
277 primdata_[p12]._0_Overlap_0_y[0] *
278 primdata_[p12]._0_Overlap_0_z[0];
280 case Operator::nuclear:
281 case Operator::erf_nuclear:
282 case Operator::erfc_nuclear:
283 for (
auto p12 = 0; p12 != primdata_[0].contrdepth; ++p12)
284 result += primdata_[p12].LIBINT_T_S_ELECPOT_S(0)[0];
287 assert(
false &&
"missing case in switch");
289 primdata_[0].targets[0] = &result;
291 const auto buildfnidx = s1.
contr[0].l * hard_lmax_ + s2.
contr[0].l;
292 assert(buildfnptrs_[buildfnidx] &&
"null build function ptr");
293 buildfnptrs_[buildfnidx](&primdata_[0]);
295 if (accumulate_ints_in_scratch) {
302 if (deriv_order_ <= 1) {
305 auto s_target = &scratch_[0];
306 for (
auto s = 0; s != ntargets; ++s, s_target += ncart12)
308 std::transform(primdata_[0].targets[s],
309 primdata_[0].targets[s] + ncart12, s_target,
310 s_target, std::plus<value_type>());
312 std::copy(primdata_[0].targets[s],
313 primdata_[0].targets[s] + ncart12, s_target);
321 if (deriv_order_ > 0) {
322 switch (deriv_order_) {
328 assert(ntargets == 6 &&
"unexpected # of targets");
329 auto dest = &scratch_[0] + (6 + pset * 3) * ncart12;
330 for (
auto s = 0; s != 3; ++s, dest += ncart12) {
331 auto src = primdata_[0].targets[s];
332 for (
auto i = 0; i != ncart12; ++i) {
337 for (
auto s = 3; s != 6; ++s, dest += ncart12) {
338 auto src = primdata_[0].targets[s];
339 for (
auto i = 0; i != ncart12; ++i) {
349 #define upper_triangle_index_ord(n2, i, j) ((i) * ((n2) - (i)-1) / 2 + (j)) 351 #define upper_triangle_index(n2, i, j) \ 352 upper_triangle_index_ord(n2, std::min((i), (j)), std::max((i), (j))) 356 const auto ncoords_times_two = nderivcoord * 2;
357 for (
auto d0 = 0, d01 = 0; d0 != 6; ++d0) {
358 for (
auto d1 = d0; d1 != 6; ++d1, ++d01) {
359 const auto d01_full =
360 upper_triangle_index_ord(ncoords_times_two, d0, d1);
361 auto tgt = &scratch_[d01_full * ncart12];
363 std::transform(primdata_[0].targets[d01],
364 primdata_[0].targets[d01] + ncart12, tgt,
365 tgt, std::plus<value_type>());
367 std::copy(primdata_[0].targets[d01],
368 primdata_[0].targets[d01] + ncart12, tgt);
377 const auto c1 = 2 + pset;
378 for (
auto c0 = 0; c0 != 2; ++c0) {
379 for (
auto xyz0 = 0; xyz0 != 3; ++xyz0) {
380 const auto coord0 = c0 * 3 + xyz0;
381 for (
auto xyz1 = 0; xyz1 != 3; ++xyz1) {
382 const auto coord1 = c1 * 3 + xyz1;
384 const auto coord01_abs = upper_triangle_index_ord(
385 ncoords_times_two, coord0, coord1);
386 auto tgt = &scratch_[coord01_abs * ncart12];
390 auto coord1_A = xyz1;
391 const auto coord01_A =
392 upper_triangle_index(12, coord0, coord1_A);
393 const auto src = primdata_[0].targets[coord01_A];
394 for (
auto i = 0; i != ncart12; ++i) tgt[i] = -src[i];
399 auto coord1_B = 3 + xyz1;
400 const auto coord01_B =
401 upper_triangle_index(12, coord0, coord1_B);
402 const auto src = primdata_[0].targets[coord01_B];
403 for (
auto i = 0; i != ncart12; ++i) tgt[i] -= src[i];
411 const auto c0 = 2 + pset;
413 for (
auto xyz0 = 0; xyz0 != 3; ++xyz0) {
414 const auto coord0 = c0 * 3 + xyz0;
415 for (
auto xyz1 = xyz0; xyz1 != 3; ++xyz1) {
416 const auto coord1 = c1 * 3 + xyz1;
418 const auto coord01_abs = upper_triangle_index_ord(
419 ncoords_times_two, coord0, coord1);
420 auto tgt = &scratch_[coord01_abs * ncart12];
424 auto coord0_A = xyz0;
425 auto coord1_A = xyz1;
426 const auto coord01_AA =
427 upper_triangle_index_ord(12, coord0_A, coord1_A);
428 const auto src = primdata_[0].targets[coord01_AA];
429 for (
auto i = 0; i != ncart12; ++i) tgt[i] = src[i];
434 auto coord0_A = xyz0;
435 auto coord1_B = 3 + xyz1;
436 const auto coord01_AB =
437 upper_triangle_index_ord(12, coord0_A, coord1_B);
438 const auto src = primdata_[0].targets[coord01_AB];
439 for (
auto i = 0; i != ncart12; ++i) tgt[i] += src[i];
444 auto coord0_B = 3 + xyz0;
445 auto coord1_A = xyz1;
446 const auto coord01_BA =
447 upper_triangle_index_ord(12, coord1_A, coord0_B);
448 const auto src = primdata_[0].targets[coord01_BA];
449 for (
auto i = 0; i != ncart12; ++i) tgt[i] += src[i];
454 auto coord0_B = 3 + xyz0;
455 auto coord1_B = 3 + xyz1;
456 const auto coord01_BB =
457 upper_triangle_index_ord(12, coord0_B, coord1_B);
458 const auto src = primdata_[0].targets[coord01_BB];
459 for (
auto i = 0; i != ncart12; ++i) tgt[i] += src[i];
465 #undef upper_triangle_index 469 assert(deriv_order_ <= 2 &&
"feature not implemented");
474 using ShellSetDerivIterator =
476 std::vector<unsigned int>>;
477 ShellSetDerivIterator shellset_gaussian_diter(deriv_order_, 2);
478 ShellSetDerivIterator shellset_full_diter(deriv_order_,
479 nderivcenters_shset);
480 std::vector<unsigned int> full_deriv(3 * nderivcenters_shset, 0);
482 while (shellset_gaussian_diter) {
484 const auto& s1s2_deriv = *shellset_gaussian_diter;
485 std::copy(std::begin(s1s2_deriv), std::end(s1s2_deriv),
486 std::begin(full_deriv));
487 const auto full_rank = ShellSetDerivIterator::rank(full_deriv);
488 targets_[full_rank] = primdata_[0].targets[s];
500 if (tform_to_solids) {
503 auto* spherical_ints =
504 (accumulate_ints_in_scratch) ? scratch2_ : &scratch_[0];
508 for (
auto s = 0ul; s != num_shellsets_computed;
509 ++s, spherical_ints += n12) {
510 auto cartesian_ints = accumulate_ints_in_scratch
511 ? &scratch_[s * ncart12]
512 : primdata_[0].targets[s];
515 libint2::solidharmonics::tform(l1, l2, cartesian_ints, spherical_ints);
517 if (s1.
contr[0].pure)
518 libint2::solidharmonics::tform_rows(l1, n2, cartesian_ints,
521 libint2::solidharmonics::tform_cols(n1, l2, cartesian_ints,
525 targets_[s] = spherical_ints;
530 for (
auto s = 0ul; s != num_shellsets_computed; ++s) {
531 auto cartesian_ints = accumulate_ints_in_scratch
532 ? &scratch_[s * ncart12]
533 : primdata_[0].targets[s];
534 targets_[s] = cartesian_ints;
538 if (cartesian_shell_normalization() == CartesianShellNormalization::uniform) {
540 for (
auto s = 0ul; s != num_shellsets_computed; ++s) {
549 __libint2_engine_inline
void Engine::_initialize() {
550 #define BOOST_PP_NBODYENGINE_MCR3_ncenter(product) \ 551 BOOST_PP_TUPLE_ELEM(3, 1, product) 553 #define BOOST_PP_NBODYENGINE_MCR3_default_ncenter(product) \ 554 BOOST_PP_IIF(BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(3, 0, product), \ 555 BOOST_PP_NBODY_OPERATOR_LAST_ONEBODY_INDEX), \ 558 #define BOOST_PP_NBODYENGINE_MCR3_NCENTER(product) \ 560 BOOST_PP_NOT_EQUAL(BOOST_PP_NBODYENGINE_MCR3_ncenter(product), \ 561 BOOST_PP_NBODYENGINE_MCR3_default_ncenter(product)), \ 562 BOOST_PP_NBODYENGINE_MCR3_ncenter(product), BOOST_PP_EMPTY()) 564 #define BOOST_PP_NBODYENGINE_MCR3_OPER(product) \ 565 BOOST_PP_LIST_AT(BOOST_PP_NBODY_OPERATOR_LIST, \ 566 BOOST_PP_TUPLE_ELEM(3, 0, product)) 568 #define BOOST_PP_NBODYENGINE_MCR3_DERIV(product) \ 569 BOOST_PP_IIF(BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(3, 2, product), 0), \ 570 BOOST_PP_TUPLE_ELEM(3, 2, product), BOOST_PP_EMPTY()) 572 #define BOOST_PP_NBODYENGINE_MCR3_task(product) \ 573 BOOST_PP_CAT(BOOST_PP_CAT(BOOST_PP_NBODYENGINE_MCR3_ncenter(product), \ 574 BOOST_PP_NBODYENGINE_MCR3_OPER(product)), \ 575 BOOST_PP_NBODYENGINE_MCR3_DERIV(product)) 577 #define BOOST_PP_NBODYENGINE_MCR3_TASK(product) \ 579 BOOST_PP_CAT(LIBINT2_TASK_EXISTS_, \ 580 BOOST_PP_NBODYENGINE_MCR3_task(product)), \ 581 BOOST_PP_CAT(BOOST_PP_CAT(BOOST_PP_NBODYENGINE_MCR3_NCENTER(product), \ 582 BOOST_PP_NBODYENGINE_MCR3_OPER(product)), \ 583 BOOST_PP_NBODYENGINE_MCR3_DERIV(product)), \ 586 #define BOOST_PP_NBODYENGINE_MCR3(r, product) \ 587 if (static_cast<int>(oper_) == BOOST_PP_TUPLE_ELEM(3, 0, product) && \ 588 static_cast<int>(rank(braket_)) == BOOST_PP_TUPLE_ELEM(3, 1, product) && \ 589 deriv_order_ == BOOST_PP_TUPLE_ELEM(3, 2, product)) { \ 590 hard_lmax_ = BOOST_PP_CAT(LIBINT2_MAX_AM_, \ 591 BOOST_PP_NBODYENGINE_MCR3_TASK(product)) + \ 593 hard_default_lmax_ = \ 594 BOOST_PP_IF(BOOST_PP_IS_1(BOOST_PP_CAT(LIBINT2_CENTER_DEPENDENT_MAX_AM_, \ 595 BOOST_PP_NBODYENGINE_MCR3_task(product))), \ 596 BOOST_PP_CAT(LIBINT2_MAX_AM_, \ 597 BOOST_PP_CAT(default, \ 598 BOOST_PP_NBODYENGINE_MCR3_DERIV(product) \ 600 ) + 1, std::numeric_limits<int>::max()); \ 602 BOOST_PP_IF(BOOST_PP_IS_1(BOOST_PP_CAT(LIBINT2_CENTER_DEPENDENT_MAX_AM_, \ 603 BOOST_PP_NBODYENGINE_MCR3_task(product))), \ 604 std::max(hard_lmax_,hard_default_lmax_), hard_lmax_); \ 605 if (lmax_ >= lmax) { \ 606 throw Engine::lmax_exceeded( \ 607 BOOST_PP_STRINGIZE(BOOST_PP_NBODYENGINE_MCR3_TASK(product)), \ 610 if (stack_size_ > 0) \ 611 libint2_cleanup_default(&primdata_[0]); \ 612 stack_size_ = LIBINT2_PREFIXED_NAME(BOOST_PP_CAT( \ 613 libint2_need_memory_, BOOST_PP_NBODYENGINE_MCR3_TASK(product)))( \ 615 LIBINT2_PREFIXED_NAME( \ 616 BOOST_PP_CAT(libint2_init_, BOOST_PP_NBODYENGINE_MCR3_TASK(product))) \ 617 (&primdata_[0], lmax_, 0); \ 618 BOOST_PP_IF(BOOST_PP_IS_1(LIBINT2_FLOP_COUNT), \ 619 LIBINT2_PREFIXED_NAME(libint2_init_flopcounter) \ 620 (&primdata_[0], primdata_.size()), BOOST_PP_EMPTY()); \ 621 buildfnptrs_ = to_ptr1(LIBINT2_PREFIXED_NAME(BOOST_PP_CAT( \ 622 libint2_build_, BOOST_PP_NBODYENGINE_MCR3_TASK(product)))); \ 627 BOOST_PP_LIST_FOR_EACH_PRODUCT(
628 BOOST_PP_NBODYENGINE_MCR3, 3,
629 (BOOST_PP_NBODY_OPERATOR_INDEX_LIST, BOOST_PP_NBODY_BRAKET_RANK_LIST,
630 BOOST_PP_NBODY_DERIV_ORDER_LIST))
634 "missing case in switch");
637 __libint2_engine_inline
void Engine::initialize(
size_t max_nprim) {
639 assert(deriv_order_ <= LIBINT2_MAX_DERIV_ORDER &&
640 "exceeded the max derivative order of the library");
643 #ifndef INCLUDE_ONEBODY 644 assert(braket_ != BraKet::x_x &&
645 "this braket type not supported by the library; give --enable-1body to configure");
648 assert(braket_ != BraKet::xx_xx &&
649 "this braket type not supported by the library; give --enable-eri to configure");
652 assert((braket_ != BraKet::xs_xx && braket_ != BraKet::xx_xs) &&
653 "this braket type not supported by the library; give --enable-eri3 to configure");
656 assert(braket_ != BraKet::xs_xs &&
657 "this braket type not supported by the library; give --enable-eri2 to configure");
662 throw using_default_initialized();
665 if (braket_ == BraKet::invalid) braket_ = default_braket(oper_);
667 if (max_nprim != 0) primdata_.resize(std::pow(max_nprim, braket_rank()));
671 decltype(targets_)::allocator_type alloc(primdata_[0].targets);
672 targets_ = decltype(targets_)(alloc);
678 const auto permutable_targets =
680 (braket_ == BraKet::xx_xx || braket_ == BraKet::xs_xx ||
681 braket_ == BraKet::xx_xs);
682 if (permutable_targets)
683 targets_.reserve(max_ntargets + 1);
685 targets_.reserve(max_ntargets);
689 #ifdef LIBINT2_ENGINE_TIMERS 690 timers.set_now_overhead(25);
692 #ifdef LIBINT2_PROFILE 693 primdata_[0].timers->set_now_overhead(25);
700 __libint2_engine_inline std::vector<Engine::compute2_ptr_type>
701 init_compute2_ptrs() {
702 auto max_ncompute2_ptrs = nopers_2body * nbrakets_2body * nderivorders_2body;
703 std::vector<Engine::compute2_ptr_type> result(max_ncompute2_ptrs,
nullptr);
705 #define BOOST_PP_NBODYENGINE_MCR7(r, product) \ 706 if (BOOST_PP_TUPLE_ELEM(3, 0, product) >= \ 707 static_cast<int>(Operator::first_2body_oper) && \ 708 BOOST_PP_TUPLE_ELEM(3, 0, product) <= \ 709 static_cast<int>(Operator::last_2body_oper) && \ 710 BOOST_PP_TUPLE_ELEM(3, 1, product) >= \ 711 static_cast<int>(BraKet::first_2body_braket) && \ 712 BOOST_PP_TUPLE_ELEM(3, 1, product) <= \ 713 static_cast<int>(BraKet::last_2body_braket)) { \ 714 auto compute_ptr_idx = ((BOOST_PP_TUPLE_ELEM(3, 0, product) - \ 715 static_cast<int>(Operator::first_2body_oper)) * \ 717 (BOOST_PP_TUPLE_ELEM(3, 1, product) - \ 718 static_cast<int>(BraKet::first_2body_braket))) * \ 719 nderivorders_2body + \ 720 BOOST_PP_TUPLE_ELEM(3, 2, product); \ 721 result.at(compute_ptr_idx) = &Engine::compute2< \ 722 static_cast<Operator>(BOOST_PP_TUPLE_ELEM(3, 0, product)), \ 723 static_cast<BraKet>(BOOST_PP_TUPLE_ELEM(3, 1, product)), \ 724 BOOST_PP_TUPLE_ELEM(3, 2, product)>; \ 727 BOOST_PP_LIST_FOR_EACH_PRODUCT(
728 BOOST_PP_NBODYENGINE_MCR7, 3,
729 (BOOST_PP_NBODY_OPERATOR_INDEX_LIST, BOOST_PP_NBODY_BRAKET_INDEX_LIST,
730 BOOST_PP_NBODY_DERIV_ORDER_LIST))
736 __libint2_engine_inline
const std::vector<Engine::compute2_ptr_type>&
737 Engine::compute2_ptrs()
const {
738 static std::vector<compute2_ptr_type> compute2_ptrs_ =
739 detail::init_compute2_ptrs();
740 return compute2_ptrs_;
743 __libint2_engine_inline
unsigned int Engine::nparams()
const {
745 case Operator::nuclear:
746 return any_cast<
const operator_traits<Operator::nuclear>::oper_params_type&>(params_)
748 case Operator::erf_nuclear:
749 case Operator::erfc_nuclear:
750 return std::get<1>(any_cast<
const operator_traits<Operator::erfc_nuclear>::oper_params_type&>(params_))
757 __libint2_engine_inline
unsigned int Engine::nopers()
const {
758 switch (static_cast<int>(oper_)) {
759 #define BOOST_PP_NBODYENGINE_MCR4(r, data, i, elem) \ 761 return operator_traits<static_cast<Operator>(i)>::nopers; 762 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR4, _,
763 BOOST_PP_NBODY_OPERATOR_LIST)
767 assert(
false &&
"missing case in switch");
772 __libint2_engine_inline any Engine::enforce_params_type<any>(
773 Operator oper,
const any& params,
bool throw_if_wrong_type) {
775 switch (static_cast<int>(oper)) {
776 #define BOOST_PP_NBODYENGINE_MCR5A(r, data, i, elem) \ 778 if (any_cast<operator_traits<static_cast<Operator>(i)>::oper_params_type>( \ 779 ¶ms) != nullptr) { \ 782 if (throw_if_wrong_type) throw bad_any_cast(); \ 783 result = operator_traits<static_cast<Operator>(i)>::default_params(); \ 787 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR5A, _,
788 BOOST_PP_NBODY_OPERATOR_LIST)
791 assert(
false &&
"missing case in switch");
796 template <
typename Params>
797 __libint2_engine_inline any Engine::enforce_params_type(
798 Operator oper,
const Params& params,
bool throw_if_wrong_type) {
800 switch (static_cast<int>(oper)) {
801 #define BOOST_PP_NBODYENGINE_MCR5B(r, data, i, elem) \ 803 if (std::is_same<Params, operator_traits<static_cast<Operator>( \ 804 i)>::oper_params_type>::value) { \ 807 if (throw_if_wrong_type) throw std::bad_cast(); \ 808 result = operator_traits<static_cast<Operator>(i)>::default_params(); \ 812 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR5B, _,
813 BOOST_PP_NBODY_OPERATOR_LIST)
816 assert(
false &&
"missing case in switch");
821 __libint2_engine_inline any Engine::make_core_eval_pack(Operator oper)
const {
823 switch (static_cast<int>(oper)) {
824 #define BOOST_PP_NBODYENGINE_MCR6(r, data, i, elem) \ 826 result = libint2::detail::make_compressed_pair( \ 827 operator_traits<static_cast<Operator>(i)>::core_eval_type::instance( \ 828 braket_rank() * lmax_ + deriv_order_, \ 829 std::numeric_limits<scalar_type>::epsilon()), \ 830 libint2::detail::CoreEvalScratch< \ 831 operator_traits<static_cast<Operator>(i)>::core_eval_type>( \ 832 braket_rank() * lmax_ + deriv_order_)); \ 833 assert(any_cast<detail::core_eval_pack_type<static_cast<Operator>(i)>>( \ 834 &result) != nullptr); \ 837 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR6, _,
838 BOOST_PP_NBODY_OPERATOR_LIST)
841 assert(
false &&
"missing case in switch");
846 __libint2_engine_inline
void Engine::init_core_ints_params(
const any& params) {
847 if (oper_ == Operator::delcgtg2) {
853 any_cast<
const operator_traits<Operator::delcgtg2>::oper_params_type&>(params);
854 const auto ng = oparams.size();
855 operator_traits<Operator::delcgtg2>::oper_params_type core_ints_params;
856 core_ints_params.reserve(ng * (ng + 1) / 2);
857 for (
size_t b = 0; b < ng; ++b)
858 for (
size_t k = 0; k <= b; ++k) {
859 const auto gexp = oparams[b].first + oparams[k].first;
860 const auto gcoeff = oparams[b].second * oparams[k].second *
862 const auto gcoeff_rescaled =
863 4 * oparams[b].first * oparams[k].first * gcoeff;
864 core_ints_params.push_back(std::make_pair(gexp, gcoeff_rescaled));
866 core_ints_params_ = core_ints_params;
868 core_ints_params_ = params;
872 __libint2_engine_inline
void Engine::compute_primdata(Libint_t& primdata,
const Shell& s1,
873 const Shell& s2,
size_t p1,
size_t p2,
875 const auto& A = s1.O;
876 const auto& B = s2.O;
878 const auto alpha1 = s1.alpha[p1];
879 const auto alpha2 = s2.alpha[p2];
881 const auto c1 = s1.contr[0].coeff[p1];
882 const auto c2 = s2.contr[0].coeff[p2];
884 const auto gammap = alpha1 + alpha2;
885 const auto oogammap = 1 / gammap;
886 const auto rhop_over_alpha1 = alpha2 * oogammap;
887 const auto rhop = alpha1 * rhop_over_alpha1;
888 const auto Px = (alpha1 * A[0] + alpha2 * B[0]) * oogammap;
889 const auto Py = (alpha1 * A[1] + alpha2 * B[1]) * oogammap;
890 const auto Pz = (alpha1 * A[2] + alpha2 * B[2]) * oogammap;
891 const auto AB_x = A[0] - B[0];
892 const auto AB_y = A[1] - B[1];
893 const auto AB_z = A[2] - B[2];
894 const auto AB2_x = AB_x * AB_x;
895 const auto AB2_y = AB_y * AB_y;
896 const auto AB2_z = AB_z * AB_z;
898 assert(LIBINT2_SHELLQUARTET_SET == LIBINT2_SHELLQUARTET_SET_STANDARD &&
"non-standard shell ordering");
900 const auto oper_is_nuclear =
901 (oper_ == Operator::nuclear || oper_ == Operator::erf_nuclear ||
902 oper_ == Operator::erfc_nuclear);
905 const auto l1 = s1.contr[0].l;
906 const auto l2 = s2.contr[0].l;
907 const bool use_hrr = (oper_is_nuclear || oper_ == Operator::sphemultipole) && l1 > 0 && l2 > 0;
909 const bool hrr_ket_to_bra = l1 >= l2;
911 if (hrr_ket_to_bra) {
912 #if LIBINT2_DEFINED(eri, AB_x) 913 primdata.AB_x[0] = AB_x;
915 #if LIBINT2_DEFINED(eri, AB_y) 916 primdata.AB_y[0] = AB_y;
918 #if LIBINT2_DEFINED(eri, AB_z) 919 primdata.AB_z[0] = AB_z;
923 #if LIBINT2_DEFINED(eri, BA_x) 924 primdata.BA_x[0] = - AB_x;
926 #if LIBINT2_DEFINED(eri, BA_y) 927 primdata.BA_y[0] = - AB_y;
929 #if LIBINT2_DEFINED(eri, BA_z) 930 primdata.BA_z[0] = - AB_z;
937 #if LIBINT2_DEFINED(eri, PA_x) 938 primdata.PA_x[0] = Px - A[0];
940 #if LIBINT2_DEFINED(eri, PA_y) 941 primdata.PA_y[0] = Py - A[1];
943 #if LIBINT2_DEFINED(eri, PA_z) 944 primdata.PA_z[0] = Pz - A[2];
949 #if LIBINT2_DEFINED(eri, PB_x) 950 primdata.PB_x[0] = Px - B[0];
952 #if LIBINT2_DEFINED(eri, PB_y) 953 primdata.PB_y[0] = Py - B[1];
955 #if LIBINT2_DEFINED(eri, PB_z) 956 primdata.PB_z[0] = Pz - B[2];
960 if (oper_ == Operator::emultipole1 || oper_ == Operator::emultipole2 ||
961 oper_ == Operator::emultipole3) {
962 const auto& O = any_cast<
const operator_traits<
963 Operator::emultipole1>::oper_params_type&>(params_);
964 #if LIBINT2_DEFINED(eri, BO_x) 965 primdata.BO_x[0] = B[0] - O[0];
967 #if LIBINT2_DEFINED(eri, BO_y) 968 primdata.BO_y[0] = B[1] - O[1];
970 #if LIBINT2_DEFINED(eri, BO_z) 971 primdata.BO_z[0] = B[2] - O[2];
974 if (oper_ == Operator::sphemultipole) {
975 const auto& O = any_cast<
const operator_traits<
976 Operator::emultipole1>::oper_params_type&>(params_);
977 #if LIBINT2_DEFINED(eri, PO_x) 978 primdata.PO_x[0] = Px - O[0];
980 #if LIBINT2_DEFINED(eri, PO_y) 981 primdata.PO_y[0] = Py - O[1];
983 #if LIBINT2_DEFINED(eri, PO_z) 984 primdata.PO_z[0] = Pz - O[2];
986 #if LIBINT2_DEFINED(eri, PO2) 987 primdata.PO2[0] = (Px - O[0])*(Px - O[0]) + (Py - O[1])*(Py - O[1]) + (Pz - O[2])*(Pz - O[2]);
991 #if LIBINT2_DEFINED(eri, oo2z) 992 primdata.oo2z[0] = 0.5 * oogammap;
995 decltype(c1) sqrt_PI(1.77245385090551602729816748334);
996 const auto xyz_pfac = sqrt_PI * sqrt(oogammap);
997 const auto ovlp_ss_x = exp(-rhop * AB2_x) * xyz_pfac * c1 * c2;
998 const auto ovlp_ss_y = exp(-rhop * AB2_y) * xyz_pfac;
999 const auto ovlp_ss_z = exp(-rhop * AB2_z) * xyz_pfac;
1001 primdata._0_Overlap_0_x[0] = ovlp_ss_x;
1002 primdata._0_Overlap_0_y[0] = ovlp_ss_y;
1003 primdata._0_Overlap_0_z[0] = ovlp_ss_z;
1005 if (oper_ == Operator::kinetic || (deriv_order_ > 0)) {
1006 #if LIBINT2_DEFINED(eri, two_alpha0_bra) 1007 primdata.two_alpha0_bra[0] = 2.0 * alpha1;
1009 #if LIBINT2_DEFINED(eri, two_alpha0_ket) 1010 primdata.two_alpha0_ket[0] = 2.0 * alpha2;
1014 if (oper_is_nuclear) {
1016 const auto& params = (oper_ == Operator::nuclear) ?
1017 any_cast<
const operator_traits<Operator::nuclear>::oper_params_type&>(params_) :
1018 std::get<1>(any_cast<
const operator_traits<Operator::erfc_nuclear>::oper_params_type&>(params_));
1020 const auto& C = params[oset].second;
1021 const auto& q = params[oset].first;
1022 #if LIBINT2_DEFINED(eri, PC_x) 1023 primdata.PC_x[0] = Px - C[0];
1025 #if LIBINT2_DEFINED(eri, PC_y) 1026 primdata.PC_y[0] = Py - C[1];
1028 #if LIBINT2_DEFINED(eri, PC_z) 1029 primdata.PC_z[0] = Pz - C[2];
1032 #if LIBINT2_DEFINED(eri, rho12_over_alpha1) || \ 1033 LIBINT2_DEFINED(eri, rho12_over_alpha2) 1034 if (deriv_order_ > 0) {
1035 #if LIBINT2_DEFINED(eri, rho12_over_alpha1) 1036 primdata.rho12_over_alpha1[0] = rhop_over_alpha1;
1038 #if LIBINT2_DEFINED(eri, rho12_over_alpha2) 1039 primdata.rho12_over_alpha2[0] = alpha1 * oogammap;
1043 #if LIBINT2_DEFINED(eri, PC_x) && LIBINT2_DEFINED(eri, PC_y) && \ 1044 LIBINT2_DEFINED(eri, PC_z) 1045 const auto PC2 = primdata.PC_x[0] * primdata.PC_x[0] +
1046 primdata.PC_y[0] * primdata.PC_y[0] +
1047 primdata.PC_z[0] * primdata.PC_z[0];
1048 const scalar_type U = gammap * PC2;
1049 const scalar_type rho = rhop;
1050 const auto mmax = s1.contr[0].l + s2.contr[0].l + deriv_order_;
1051 auto* fm_ptr = &(primdata.LIBINT_T_S_ELECPOT_S(0)[0]);
1052 if (oper_ == Operator::nuclear) {
1053 auto fm_engine_ptr =
1054 any_cast<
const detail::core_eval_pack_type<Operator::nuclear>&>(core_eval_pack_)
1056 fm_engine_ptr->eval(fm_ptr, U, mmax);
1057 }
else if (oper_ == Operator::erf_nuclear) {
1058 const auto& core_eval_ptr =
1059 any_cast<
const detail::core_eval_pack_type<Operator::erf_nuclear>&>(core_eval_pack_)
1061 auto core_ints_params =
1062 std::get<0>(any_cast<
const typename operator_traits<
1063 Operator::erf_nuclear>::oper_params_type&>(core_ints_params_));
1064 core_eval_ptr->eval(fm_ptr, rho, U, mmax, core_ints_params);
1065 }
else if (oper_ == Operator::erfc_nuclear) {
1066 const auto& core_eval_ptr =
1067 any_cast<
const detail::core_eval_pack_type<Operator::erfc_nuclear>&>(core_eval_pack_)
1069 auto core_ints_params =
1070 std::get<0>(any_cast<
const typename operator_traits<
1071 Operator::erfc_nuclear>::oper_params_type&>(core_ints_params_));
1072 core_eval_ptr->eval(fm_ptr, rho, U, mmax, core_ints_params);
1075 decltype(U) two_o_sqrt_PI(1.12837916709551257389615890312);
1077 -q * sqrt(gammap) * two_o_sqrt_PI * ovlp_ss_x * ovlp_ss_y * ovlp_ss_z;
1078 const auto m_fence = mmax + 1;
1079 for (
auto m = 0; m != m_fence; ++m) {
1090 template <Operator op, BraKet bk,
size_t deriv_order>
1091 __libint2_engine_inline
const Engine::target_ptr_vec& Engine::compute2(
1094 const ShellPair* tspbra,
const ShellPair* tspket) {
1095 assert(op == oper_ &&
"Engine::compute2 -- operator mismatch");
1096 assert(bk == braket_ &&
"Engine::compute2 -- braket mismatch");
1097 assert(deriv_order == deriv_order_ &&
1098 "Engine::compute2 -- deriv_order mismatch");
1099 assert(((tspbra ==
nullptr && tspket ==
nullptr) || (tspbra !=
nullptr && tspket !=
nullptr)) &&
1100 "Engine::compute2 -- expects zero or two ShellPair objects");
1107 assert((tbra1.ncontr() == 1 && tbra2.ncontr() == 1 && tket1.ncontr() == 1 &&
1108 tket2.ncontr() == 1) &&
"generally-contracted shells are not yet supported");
1111 assert(tbra1.
contr[0].l <= lmax_ &&
"the angular momentum limit is exceeded");
1112 assert(tbra2.
contr[0].l <= lmax_ &&
"the angular momentum limit is exceeded");
1113 assert(tket1.
contr[0].l <= lmax_ &&
"the angular momentum limit is exceeded");
1114 assert(tket2.
contr[0].l <= lmax_ &&
"the angular momentum limit is exceeded");
1116 #if LIBINT2_SHELLQUARTET_SET == \ 1117 LIBINT2_SHELLQUARTET_SET_STANDARD // standard angular momentum ordering 1118 const auto swap_tbra = (tbra1.
contr[0].l < tbra2.
contr[0].l);
1119 const auto swap_tket = (tket1.
contr[0].l < tket2.
contr[0].l);
1120 const auto swap_braket =
1121 ((braket_ == BraKet::xx_xx) && (tbra1.
contr[0].l + tbra2.
contr[0].l >
1123 braket_ == BraKet::xx_xs;
1124 #else // orca angular momentum ordering 1125 const auto swap_tbra = (tbra1.
contr[0].l > tbra2.
contr[0].l);
1126 const auto swap_tket = (tket1.
contr[0].l > tket2.
contr[0].l);
1127 const auto swap_braket =
1128 ((braket_ == BraKet::xx_xx) && (tbra1.
contr[0].l + tbra2.
contr[0].l <
1130 braket_ == BraKet::xx_xs;
1131 assert(
false &&
"feature not implemented");
1134 swap_braket ? (swap_tket ? tket2 : tket1) : (swap_tbra ? tbra2 : tbra1);
1136 swap_braket ? (swap_tket ? tket1 : tket2) : (swap_tbra ? tbra1 : tbra2);
1138 swap_braket ? (swap_tbra ? tbra2 : tbra1) : (swap_tket ? tket2 : tket1);
1140 swap_braket ? (swap_tbra ? tbra1 : tbra2) : (swap_tket ? tket1 : tket2);
1141 const auto swap_bra = swap_braket ? swap_tket : swap_tbra;
1142 const auto swap_ket = swap_braket ? swap_tbra : swap_tket;
1144 const auto* spbra_precomputed = swap_braket ? tspket : tspbra;
1145 const auto* spket_precomputed = swap_braket ? tspbra : tspket;
1147 const auto tform = bra1.contr[0].pure || bra2.contr[0].pure ||
1148 ket1.contr[0].pure || ket2.contr[0].pure;
1149 const auto permute = swap_braket || swap_tbra || swap_tket;
1150 const auto use_scratch = permute || tform;
1153 auto nprim_bra1 = bra1.nprim();
1154 auto nprim_bra2 = bra2.nprim();
1155 auto nprim_ket1 = ket1.nprim();
1156 auto nprim_ket2 = ket2.nprim();
1159 auto lmax = std::max(std::max(bra1.contr[0].l, bra2.contr[0].l),
1160 std::max(ket1.contr[0].l, ket2.contr[0].l));
1161 assert(lmax <= lmax_ &&
"the angular momentum limit is exceeded");
1162 const auto lmax_bra = std::max(bra1.contr[0].l, bra2.contr[0].l);
1163 const auto lmax_ket = std::max(ket1.contr[0].l, ket2.contr[0].l);
1165 #ifdef LIBINT2_ENGINE_PROFILE_CLASS 1166 class_id id(bra1.contr[0].l, bra2.contr[0].l, ket1.contr[0].l,
1168 if (class_profiles.find(
id) == class_profiles.end()) {
1169 class_profile dummy;
1170 class_profiles[id] = dummy;
1175 #ifdef LIBINT2_ENGINE_TIMERS 1184 const ShellPair& spbra = spbra_precomputed ? *spbra_precomputed : (spbra_.init(bra1, bra2, ln_precision_), spbra_) ;
1185 const ShellPair& spket = spket_precomputed ? *spket_precomputed : (spket_.init(ket1, ket2, ln_precision_), spket_);
1188 const auto spbra_is_swapped = spbra_precomputed ? swap_bra :
false;
1189 const auto spket_is_swapped = spket_precomputed ? swap_ket :
false;
1191 using real_t = Shell::real_t;
1194 if (spbra_is_swapped) {
1195 for(
auto xyz=0; xyz!=3; ++xyz)
1196 BA[xyz] = - spbra_precomputed->AB[xyz];
1198 const auto& AB = spbra_is_swapped ? BA : spbra.AB;
1201 if (spket_is_swapped) {
1202 for(
auto xyz=0; xyz!=3; ++xyz)
1203 DC[xyz] = - spket_precomputed->AB[xyz];
1205 const auto& CD = spket_is_swapped ? DC : spket.AB;
1207 const auto& A = bra1.O;
1208 const auto& B = bra2.O;
1209 const auto& C = ket1.O;
1210 const auto& D = ket2.O;
1213 const auto npbra = spbra.primpairs.size();
1214 const auto npket = spket.primpairs.size();
1215 for (
auto pb = 0; pb != npbra; ++pb) {
1216 for (
auto pk = 0; pk != npket; ++pk) {
1218 if (spbra.primpairs[pb].scr + spket.primpairs[pk].scr >
1220 Libint_t& primdata = primdata_[p];
1221 const auto& sbra1 = bra1;
1222 const auto& sbra2 = bra2;
1223 const auto& sket1 = ket1;
1224 const auto& sket2 = ket2;
1228 const auto& spbrapp = spbra.primpairs[pbra];
1229 const auto& spketpp = spket.primpairs[pket];
1231 const auto& pbra1 = spbra_is_swapped ? spbrapp.p2 : spbrapp.p1;
1232 const auto& pbra2 = spbra_is_swapped ? spbrapp.p1 : spbrapp.p2;
1233 const auto& pket1 = spket_is_swapped ? spketpp.p2 : spketpp.p1;
1234 const auto& pket2 = spket_is_swapped ? spketpp.p1 : spketpp.p2;
1236 const auto alpha0 = sbra1.alpha[pbra1];
1237 const auto alpha1 = sbra2.alpha[pbra2];
1238 const auto alpha2 = sket1.alpha[pket1];
1239 const auto alpha3 = sket2.alpha[pket2];
1241 const auto c0 = sbra1.contr[0].coeff[pbra1];
1242 const auto c1 = sbra2.contr[0].coeff[pbra2];
1243 const auto c2 = sket1.contr[0].coeff[pket1];
1244 const auto c3 = sket2.contr[0].coeff[pket2];
1246 const auto amtot = sbra1.contr[0].l + sket1.contr[0].l +
1247 sbra2.contr[0].l + sket2.contr[0].l;
1249 const auto gammap = alpha0 + alpha1;
1250 const auto oogammap = spbrapp.one_over_gamma;
1251 const auto rhop = alpha0 * alpha1 * oogammap;
1253 const auto gammaq = alpha2 + alpha3;
1254 const auto oogammaq = spketpp.one_over_gamma;
1255 const auto rhoq = alpha2 * alpha3 * oogammaq;
1257 const auto& P = spbrapp.P;
1258 const auto& Q = spketpp.P;
1259 const auto PQx = P[0] - Q[0];
1260 const auto PQy = P[1] - Q[1];
1261 const auto PQz = P[2] - Q[2];
1262 const auto PQ2 = PQx * PQx + PQy * PQy + PQz * PQz;
1264 const auto K12 = spbrapp.K * spketpp.K;
1265 decltype(K12) two_times_M_PI_to_25(
1266 34.986836655249725693);
1267 const auto gammapq = gammap + gammaq;
1268 const auto sqrt_gammapq = sqrt(gammapq);
1269 const auto oogammapq = 1.0 / (gammapq);
1270 auto pfac = two_times_M_PI_to_25 * K12 * sqrt_gammapq * oogammapq;
1271 pfac *= c0 * c1 * c2 * c3;
1273 if (std::abs(pfac) >= precision_) {
1274 const scalar_type rho = gammap * gammaq * oogammapq;
1275 const scalar_type T = PQ2 * rho;
1276 auto* gm_ptr = &(primdata.LIBINT_T_SS_EREP_SS(0)[0]);
1277 const auto mmax = amtot + deriv_order;
1279 if (!skip_core_ints) {
1281 case Operator::coulomb: {
1282 const auto& core_eval_ptr =
1283 any_cast<
const detail::core_eval_pack_type<Operator::coulomb>&>(core_eval_pack_)
1285 core_eval_ptr->eval(gm_ptr, T, mmax);
1287 case Operator::cgtg_x_coulomb: {
1288 const auto& core_eval_ptr =
1289 any_cast<
const detail::core_eval_pack_type<
1290 Operator::cgtg_x_coulomb>&>(core_eval_pack_)
1292 auto& core_eval_scratch = any_cast<detail::core_eval_pack_type<
1293 Operator::cgtg_x_coulomb>&>(core_eval_pack_)
1295 const auto& core_ints_params =
1296 any_cast<
const typename operator_traits<
1297 Operator::cgtg>::oper_params_type&>(core_ints_params_);
1298 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params,
1299 &core_eval_scratch);
1301 case Operator::cgtg: {
1302 const auto& core_eval_ptr =
1303 any_cast<
const detail::core_eval_pack_type<Operator::cgtg>&>(core_eval_pack_)
1305 const auto& core_ints_params =
1306 any_cast<
const typename operator_traits<
1307 Operator::cgtg>::oper_params_type&>(core_ints_params_);
1308 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params);
1310 case Operator::delcgtg2: {
1311 const auto& core_eval_ptr =
1312 any_cast<
const detail::core_eval_pack_type<Operator::delcgtg2>&>(core_eval_pack_)
1314 const auto& core_ints_params =
1315 any_cast<
const typename operator_traits<
1316 Operator::cgtg>::oper_params_type&>(core_ints_params_);
1317 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params);
1319 case Operator::delta: {
1320 const auto& core_eval_ptr =
1321 any_cast<
const detail::core_eval_pack_type<Operator::delta>&>(core_eval_pack_)
1323 core_eval_ptr->eval(gm_ptr, rho, T, mmax);
1325 case Operator::r12: {
1326 const auto& core_eval_ptr =
1327 any_cast<
const detail::core_eval_pack_type<Operator::r12>&>(core_eval_pack_)
1329 core_eval_ptr->eval(gm_ptr, rho, T, mmax);
1331 case Operator::erf_coulomb: {
1332 const auto& core_eval_ptr =
1333 any_cast<
const detail::core_eval_pack_type<Operator::erf_coulomb>&>(core_eval_pack_)
1335 auto core_ints_params =
1336 any_cast<
const typename operator_traits<
1337 Operator::erf_coulomb>::oper_params_type&>(core_ints_params_);
1338 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params);
1340 case Operator::erfc_coulomb: {
1341 const auto& core_eval_ptr =
1342 any_cast<
const detail::core_eval_pack_type<Operator::erfc_coulomb>&>(core_eval_pack_)
1344 auto core_ints_params =
1345 any_cast<
const typename operator_traits<
1346 Operator::erfc_coulomb>::oper_params_type&>(core_ints_params_);
1347 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params);
1349 case Operator::stg: {
1350 const auto& core_eval_ptr =
1351 any_cast<
const detail::core_eval_pack_type<Operator::stg>&>(core_eval_pack_)
1354 any_cast<
const typename operator_traits<
1355 Operator::stg>::oper_params_type&>(core_ints_params_);
1356 const auto one_over_rho = gammapq * oogammap * oogammaq;
1357 core_eval_ptr->eval_slater(gm_ptr, one_over_rho, T, mmax, zeta);
1359 case Operator::stg_x_coulomb: {
1360 const auto& core_eval_ptr =
1361 any_cast<
const detail::core_eval_pack_type<Operator::stg>&>(core_eval_pack_)
1364 any_cast<
const typename operator_traits<
1365 Operator::stg>::oper_params_type&>(core_ints_params_);
1366 const auto one_over_rho = gammapq * oogammap * oogammaq;
1367 core_eval_ptr->eval_yukawa(gm_ptr, one_over_rho, T, mmax, zeta);
1370 assert(
false &&
"missing case in a switch");
1374 for (
auto m = 0; m != mmax + 1; ++m) {
1379 if (braket_ == BraKet::xx_xx) {
1380 #if LIBINT2_DEFINED(eri, PA_x) 1381 primdata.PA_x[0] = P[0] - A[0];
1383 #if LIBINT2_DEFINED(eri, PA_y) 1384 primdata.PA_y[0] = P[1] - A[1];
1386 #if LIBINT2_DEFINED(eri, PA_z) 1387 primdata.PA_z[0] = P[2] - A[2];
1389 #if LIBINT2_DEFINED(eri, PB_x) 1390 primdata.PB_x[0] = P[0] - B[0];
1392 #if LIBINT2_DEFINED(eri, PB_y) 1393 primdata.PB_y[0] = P[1] - B[1];
1395 #if LIBINT2_DEFINED(eri, PB_z) 1396 primdata.PB_z[0] = P[2] - B[2];
1400 if (braket_ != BraKet::xs_xs) {
1401 #if LIBINT2_DEFINED(eri, QC_x) 1402 primdata.QC_x[0] = Q[0] - C[0];
1404 #if LIBINT2_DEFINED(eri, QC_y) 1405 primdata.QC_y[0] = Q[1] - C[1];
1407 #if LIBINT2_DEFINED(eri, QC_z) 1408 primdata.QC_z[0] = Q[2] - C[2];
1410 #if LIBINT2_DEFINED(eri, QD_x) 1411 primdata.QD_x[0] = Q[0] - D[0];
1413 #if LIBINT2_DEFINED(eri, QD_y) 1414 primdata.QD_y[0] = Q[1] - D[1];
1416 #if LIBINT2_DEFINED(eri, QD_z) 1417 primdata.QD_z[0] = Q[2] - D[2];
1421 if (braket_ == BraKet::xx_xx) {
1422 #if LIBINT2_DEFINED(eri, AB_x) 1423 primdata.AB_x[0] = AB[0];
1425 #if LIBINT2_DEFINED(eri, AB_y) 1426 primdata.AB_y[0] = AB[1];
1428 #if LIBINT2_DEFINED(eri, AB_z) 1429 primdata.AB_z[0] = AB[2];
1431 #if LIBINT2_DEFINED(eri, BA_x) 1432 primdata.BA_x[0] = -AB[0];
1434 #if LIBINT2_DEFINED(eri, BA_y) 1435 primdata.BA_y[0] = -AB[1];
1437 #if LIBINT2_DEFINED(eri, BA_z) 1438 primdata.BA_z[0] = -AB[2];
1442 if (braket_ != BraKet::xs_xs) {
1443 #if LIBINT2_DEFINED(eri, CD_x) 1444 primdata.CD_x[0] = CD[0];
1446 #if LIBINT2_DEFINED(eri, CD_y) 1447 primdata.CD_y[0] = CD[1];
1449 #if LIBINT2_DEFINED(eri, CD_z) 1450 primdata.CD_z[0] = CD[2];
1452 #if LIBINT2_DEFINED(eri, DC_x) 1453 primdata.DC_x[0] = -CD[0];
1455 #if LIBINT2_DEFINED(eri, DC_y) 1456 primdata.DC_y[0] = -CD[1];
1458 #if LIBINT2_DEFINED(eri, DC_z) 1459 primdata.DC_z[0] = -CD[2];
1463 const auto gammap_o_gammapgammaq = oogammapq * gammap;
1464 const auto gammaq_o_gammapgammaq = oogammapq * gammaq;
1467 (gammap_o_gammapgammaq * P[0] + gammaq_o_gammapgammaq * Q[0]);
1469 (gammap_o_gammapgammaq * P[1] + gammaq_o_gammapgammaq * Q[1]);
1471 (gammap_o_gammapgammaq * P[2] + gammaq_o_gammapgammaq * Q[2]);
1473 if (deriv_order > 0 || lmax_bra > 0) {
1474 #if LIBINT2_DEFINED(eri, WP_x) 1475 primdata.WP_x[0] = Wx - P[0];
1477 #if LIBINT2_DEFINED(eri, WP_y) 1478 primdata.WP_y[0] = Wy - P[1];
1480 #if LIBINT2_DEFINED(eri, WP_z) 1481 primdata.WP_z[0] = Wz - P[2];
1484 if (deriv_order > 0 || lmax_ket > 0) {
1485 #if LIBINT2_DEFINED(eri, WQ_x) 1486 primdata.WQ_x[0] = Wx - Q[0];
1488 #if LIBINT2_DEFINED(eri, WQ_y) 1489 primdata.WQ_y[0] = Wy - Q[1];
1491 #if LIBINT2_DEFINED(eri, WQ_z) 1492 primdata.WQ_z[0] = Wz - Q[2];
1495 #if LIBINT2_DEFINED(eri, oo2z) 1496 primdata.oo2z[0] = 0.5 * oogammap;
1498 #if LIBINT2_DEFINED(eri, oo2e) 1499 primdata.oo2e[0] = 0.5 * oogammaq;
1501 #if LIBINT2_DEFINED(eri, oo2ze) 1502 primdata.oo2ze[0] = 0.5 * oogammapq;
1504 #if LIBINT2_DEFINED(eri, roz) 1505 primdata.roz[0] = rho * oogammap;
1507 #if LIBINT2_DEFINED(eri, roe) 1508 primdata.roe[0] = rho * oogammaq;
1512 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_0_x) 1513 primdata.TwoPRepITR_pfac0_0_0_x[0] =
1514 -(alpha1 * AB[0] + alpha3 * CD[0]) * oogammap;
1516 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_0_y) 1517 primdata.TwoPRepITR_pfac0_0_0_y[0] =
1518 -(alpha1 * AB[1] + alpha3 * CD[1]) * oogammap;
1520 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_0_z) 1521 primdata.TwoPRepITR_pfac0_0_0_z[0] =
1522 -(alpha1 * AB[2] + alpha3 * CD[2]) * oogammap;
1524 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_0_x) 1525 primdata.TwoPRepITR_pfac0_1_0_x[0] =
1526 -(alpha1 * AB[0] + alpha3 * CD[0]) * oogammaq;
1528 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_0_y) 1529 primdata.TwoPRepITR_pfac0_1_0_y[0] =
1530 -(alpha1 * AB[1] + alpha3 * CD[1]) * oogammaq;
1532 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_0_z) 1533 primdata.TwoPRepITR_pfac0_1_0_z[0] =
1534 -(alpha1 * AB[2] + alpha3 * CD[2]) * oogammaq;
1536 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_1_x) 1537 primdata.TwoPRepITR_pfac0_0_1_x[0] =
1538 (alpha0 * AB[0] + alpha2 * CD[0]) * oogammap;
1540 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_1_y) 1541 primdata.TwoPRepITR_pfac0_0_1_y[0] =
1542 (alpha0 * AB[1] + alpha2 * CD[1]) * oogammap;
1544 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_1_z) 1545 primdata.TwoPRepITR_pfac0_0_1_z[0] =
1546 (alpha0 * AB[2] + alpha2 * CD[2]) * oogammap;
1548 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_1_x) 1549 primdata.TwoPRepITR_pfac0_1_1_x[0] =
1550 (alpha0 * AB[0] + alpha2 * CD[0]) * oogammaq;
1552 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_1_y) 1553 primdata.TwoPRepITR_pfac0_1_1_y[0] =
1554 (alpha0 * AB[1] + alpha2 * CD[1]) * oogammaq;
1556 #if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_1_z) 1557 primdata.TwoPRepITR_pfac0_1_1_z[0] =
1558 (alpha0 * AB[2] + alpha2 * CD[2]) * oogammaq;
1560 #if LIBINT2_DEFINED(eri, eoz) 1561 primdata.eoz[0] = gammaq * oogammap;
1563 #if LIBINT2_DEFINED(eri, zoe) 1564 primdata.zoe[0] = gammap * oogammaq;
1568 if (deriv_order > 0) {
1569 #if LIBINT2_DEFINED(eri, alpha1_rho_over_zeta2) 1570 primdata.alpha1_rho_over_zeta2[0] =
1571 alpha0 * (oogammap * gammaq_o_gammapgammaq);
1573 #if LIBINT2_DEFINED(eri, alpha2_rho_over_zeta2) 1574 primdata.alpha2_rho_over_zeta2[0] =
1575 alpha1 * (oogammap * gammaq_o_gammapgammaq);
1577 #if LIBINT2_DEFINED(eri, alpha3_rho_over_eta2) 1578 primdata.alpha3_rho_over_eta2[0] =
1579 alpha2 * (oogammaq * gammap_o_gammapgammaq);
1581 #if LIBINT2_DEFINED(eri, alpha4_rho_over_eta2) 1582 primdata.alpha4_rho_over_eta2[0] =
1583 alpha3 * (oogammaq * gammap_o_gammapgammaq);
1585 #if LIBINT2_DEFINED(eri, alpha1_over_zetapluseta) 1586 primdata.alpha1_over_zetapluseta[0] = alpha0 * oogammapq;
1588 #if LIBINT2_DEFINED(eri, alpha2_over_zetapluseta) 1589 primdata.alpha2_over_zetapluseta[0] = alpha1 * oogammapq;
1591 #if LIBINT2_DEFINED(eri, alpha3_over_zetapluseta) 1592 primdata.alpha3_over_zetapluseta[0] = alpha2 * oogammapq;
1594 #if LIBINT2_DEFINED(eri, alpha4_over_zetapluseta) 1595 primdata.alpha4_over_zetapluseta[0] = alpha3 * oogammapq;
1597 #if LIBINT2_DEFINED(eri, rho12_over_alpha1) 1598 primdata.rho12_over_alpha1[0] = alpha1 * oogammap;
1600 #if LIBINT2_DEFINED(eri, rho12_over_alpha2) 1601 primdata.rho12_over_alpha2[0] = alpha0 * oogammap;
1603 #if LIBINT2_DEFINED(eri, rho34_over_alpha3) 1604 primdata.rho34_over_alpha3[0] = alpha3 * oogammaq;
1606 #if LIBINT2_DEFINED(eri, rho34_over_alpha4) 1607 primdata.rho34_over_alpha4[0] = alpha2 * oogammaq;
1609 #if LIBINT2_DEFINED(eri, two_alpha0_bra) 1610 primdata.two_alpha0_bra[0] = 2.0 * alpha0;
1612 #if LIBINT2_DEFINED(eri, two_alpha0_ket) 1613 primdata.two_alpha0_ket[0] = 2.0 * alpha1;
1615 #if LIBINT2_DEFINED(eri, two_alpha1_bra) 1616 primdata.two_alpha1_bra[0] = 2.0 * alpha2;
1618 #if LIBINT2_DEFINED(eri, two_alpha1_ket) 1619 primdata.two_alpha1_ket[0] = 2.0 * alpha3;
1630 primdata_[0].contrdepth = p;
1633 #ifdef LIBINT2_ENGINE_TIMERS 1634 const auto t0 = timers.stop(0);
1635 #ifdef LIBINT2_ENGINE_PROFILE_CLASS 1636 class_profiles[id].prereqs += t0.count();
1637 if (primdata_[0].contrdepth != 0) {
1638 class_profiles[id].nshellset += 1;
1639 class_profiles[id].nprimset += primdata_[0].contrdepth;
1645 if (primdata_[0].contrdepth == 0) {
1646 targets_[0] =
nullptr;
1651 const auto compute_directly = lmax == 0 && deriv_order == 0;
1653 if (compute_directly) {
1654 #ifdef LIBINT2_ENGINE_TIMERS 1657 auto& stack = primdata_[0].stack[0];
1659 for (
auto p = 0; p != primdata_[0].contrdepth; ++p)
1660 stack += primdata_[p].LIBINT_T_SS_EREP_SS(0)[0];
1661 primdata_[0].targets[0] = primdata_[0].stack;
1662 #ifdef LIBINT2_ENGINE_TIMERS 1663 const auto t1 = timers.stop(1);
1664 #ifdef LIBINT2_ENGINE_PROFILE_CLASS 1665 class_profiles[id].build_vrr += t1.count();
1670 #ifdef LIBINT2_ENGINE_TIMERS 1671 #ifdef LIBINT2_PROFILE 1672 const auto t1_hrr_start = primdata_[0].timers->read(0);
1673 const auto t1_vrr_start = primdata_[0].timers->read(1);
1682 ((bra1.contr[0].l * hard_lmax_ + bra2.contr[0].l) * hard_lmax_ +
1688 case BraKet::xx_xs: assert(
false &&
"this braket is not supported");
break;
1689 case BraKet::xs_xx: {
1691 int ket_lmax = hard_lmax_;
1692 switch (deriv_order_) {
1694 #ifdef LIBINT2_CENTER_DEPENDENT_MAX_AM_3eri 1695 ket_lmax = hard_default_lmax_;
1699 #ifdef LIBINT2_CENTER_DEPENDENT_MAX_AM_3eri1 1700 ket_lmax = hard_default_lmax_;
1704 #ifdef LIBINT2_CENTER_DEPENDENT_MAX_AM_3eri2 1705 ket_lmax = hard_default_lmax_;
1708 default:assert(
false &&
"deriv_order>2 not yet supported");
1711 (bra1.contr[0].l * ket_lmax + ket1.contr[0].l) * ket_lmax +
1714 if (bra1.contr[0].l > 1)
1715 assert(bra1.contr[0].pure &&
1716 "library assumes a solid harmonics shell in bra of a 3-center " 1717 "2-body int, but a cartesian shell given");
1722 buildfnidx = bra1.contr[0].l * hard_lmax_ + ket1.contr[0].l;
1724 if (bra1.contr[0].l > 1)
1725 assert(bra1.contr[0].pure &&
1726 "library assumes solid harmonics shells in a 2-center " 1727 "2-body int, but a cartesian shell given in bra");
1728 if (ket1.contr[0].l > 1)
1729 assert(ket1.contr[0].pure &&
1730 "library assumes solid harmonics shells in a 2-center " 1731 "2-body int, but a cartesian shell given in bra");
1736 assert(
false &&
"invalid braket");
1739 assert(buildfnptrs_[buildfnidx] &&
"null build function ptr");
1740 buildfnptrs_[buildfnidx](&primdata_[0]);
1742 #ifdef LIBINT2_ENGINE_TIMERS 1743 const auto t1 = timers.stop(1);
1744 #ifdef LIBINT2_ENGINE_PROFILE_CLASS 1745 #ifndef LIBINT2_PROFILE 1746 class_profiles[id].build_vrr += t1.count();
1748 class_profiles[id].build_hrr += primdata_[0].timers->read(0) - t1_hrr_start;
1749 class_profiles[id].build_vrr += primdata_[0].timers->read(1) - t1_vrr_start;
1754 #ifdef LIBINT2_ENGINE_TIMERS 1758 const auto ntargets = nshellsets();
1762 constexpr
auto using_scalar_real =
sizeof(value_type) ==
sizeof(scalar_type);
1763 static_assert(using_scalar_real,
1764 "Libint2 C++11 API only supports scalar real types");
1765 typedef Eigen::Matrix<scalar_type, Eigen::Dynamic, Eigen::Dynamic,
1770 const auto nr1_cart = bra1.cartesian_size();
1771 const auto nr2_cart = bra2.cartesian_size();
1772 const auto nc1_cart = ket1.cartesian_size();
1773 const auto nc2_cart = ket2.cartesian_size();
1774 const auto ncol_cart = nc1_cart * nc2_cart;
1775 const auto n1234_cart = nr1_cart * nr2_cart * ncol_cart;
1776 const auto nr1 = bra1.size();
1777 const auto nr2 = bra2.size();
1778 const auto nc1 = ket1.size();
1779 const auto nc2 = ket2.size();
1780 const auto nrow = nr1 * nr2;
1781 const auto ncol = nc1 * nc2;
1784 const auto nr1_tgt = tbra1.size();
1785 const auto nr2_tgt = tbra2.size();
1786 const auto nc1_tgt = tket1.size();
1787 const auto nc2_tgt = tket2.size();
1788 const auto ncol_tgt = nc1_tgt * nc2_tgt;
1789 const auto n_tgt = nr1_tgt * nr2_tgt * ncol_tgt;
1791 auto hotscr = &scratch_[0];
1794 for (
auto s = 0; s != ntargets; ++s) {
1801 primdata_[0].targets[s];
1802 auto target = hotscr;
1804 if (bra1.contr[0].pure) {
1805 libint2::solidharmonics::transform_first(
1806 bra1.contr[0].l, nr2_cart * ncol_cart, source, target);
1807 std::swap(source, target);
1809 if (bra2.contr[0].pure) {
1810 libint2::solidharmonics::transform_inner(bra1.size(), bra2.contr[0].l,
1811 ncol_cart, source, target);
1812 std::swap(source, target);
1814 if (ket1.contr[0].pure) {
1815 libint2::solidharmonics::transform_inner(nrow, ket1.contr[0].l,
1816 nc2_cart, source, target);
1817 std::swap(source, target);
1819 if (ket2.contr[0].pure) {
1820 libint2::solidharmonics::transform_last(
1821 bra1.size() * bra2.size() * ket1.size(), ket2.contr[0].l, source,
1823 std::swap(source, target);
1829 const auto* src_row_ptr = source;
1830 auto tgt_ptr = target;
1833 switch (deriv_order) {
1839 case BraKet::xx_xx: {
1840 const unsigned mapDerivIndex1_xxxx[2][2][2][12] = {
1841 {{{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11},
1842 {0, 1, 2, 3, 4, 5, 9, 10, 11, 6, 7, 8}},
1843 {{3, 4, 5, 0, 1, 2, 6, 7, 8, 9, 10, 11},
1844 {3, 4, 5, 0, 1, 2, 9, 10, 11, 6, 7, 8}}},
1845 {{{6, 7, 8, 9, 10, 11, 0, 1, 2, 3, 4, 5},
1846 {9, 10, 11, 6, 7, 8, 0, 1, 2, 3, 4, 5}},
1847 {{6, 7, 8, 9, 10, 11, 3, 4, 5, 0, 1, 2},
1848 {9, 10, 11, 6, 7, 8, 3, 4, 5, 0, 1, 2}}}};
1849 s_target = mapDerivIndex1_xxxx[swap_braket][swap_tbra][swap_tket][s];
1853 case BraKet::xs_xx: {
1854 assert(swap_bra ==
false);
1855 assert(swap_braket ==
false);
1856 const unsigned mapDerivIndex1_xsxx[2][9] = {
1857 {0,1,2,3,4,5,6,7,8},
1860 s_target = mapDerivIndex1_xsxx[swap_tket][s];
1864 case BraKet::xs_xs: {
1865 assert(swap_bra ==
false);
1866 assert(swap_ket ==
false);
1867 assert(swap_braket ==
false);
1873 assert(
false &&
"this backet type not yet supported for 1st geometric derivatives");
1879 case BraKet::xx_xx: {
1880 const unsigned mapDerivIndex2_xxxx[2][2][2][78] = {
1881 {{{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
1882 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
1883 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
1884 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
1885 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
1886 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77},
1887 {0, 1, 2, 3, 4, 5, 9, 10, 11, 6, 7, 8, 12,
1888 13, 14, 15, 16, 20, 21, 22, 17, 18, 19, 23, 24, 25,
1889 26, 30, 31, 32, 27, 28, 29, 33, 34, 35, 39, 40, 41,
1890 36, 37, 38, 42, 43, 47, 48, 49, 44, 45, 46, 50, 54,
1891 55, 56, 51, 52, 53, 72, 73, 74, 60, 65, 69, 75, 76,
1892 61, 66, 70, 77, 62, 67, 71, 57, 58, 59, 63, 64, 68}},
1893 {{33, 34, 35, 3, 14, 24, 36, 37, 38, 39, 40, 41, 42,
1894 43, 4, 15, 25, 44, 45, 46, 47, 48, 49, 50, 5, 16,
1895 26, 51, 52, 53, 54, 55, 56, 0, 1, 2, 6, 7, 8,
1896 9, 10, 11, 12, 13, 17, 18, 19, 20, 21, 22, 23, 27,
1897 28, 29, 30, 31, 32, 57, 58, 59, 60, 61, 62, 63, 64,
1898 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77},
1899 {33, 34, 35, 3, 14, 24, 39, 40, 41, 36, 37, 38, 42,
1900 43, 4, 15, 25, 47, 48, 49, 44, 45, 46, 50, 5, 16,
1901 26, 54, 55, 56, 51, 52, 53, 0, 1, 2, 9, 10, 11,
1902 6, 7, 8, 12, 13, 20, 21, 22, 17, 18, 19, 23, 30,
1903 31, 32, 27, 28, 29, 72, 73, 74, 60, 65, 69, 75, 76,
1904 61, 66, 70, 77, 62, 67, 71, 57, 58, 59, 63, 64, 68}}},
1905 {{{57, 58, 59, 60, 61, 62, 6, 17, 27, 36, 44, 51, 63,
1906 64, 65, 66, 67, 7, 18, 28, 37, 45, 52, 68, 69, 70,
1907 71, 8, 19, 29, 38, 46, 53, 72, 73, 74, 9, 20, 30,
1908 39, 47, 54, 75, 76, 10, 21, 31, 40, 48, 55, 77, 11,
1909 22, 32, 41, 49, 56, 0, 1, 2, 3, 4, 5, 12, 13,
1910 14, 15, 16, 23, 24, 25, 26, 33, 34, 35, 42, 43, 50},
1911 {72, 73, 74, 60, 65, 69, 9, 20, 30, 39, 47, 54, 75,
1912 76, 61, 66, 70, 10, 21, 31, 40, 48, 55, 77, 62, 67,
1913 71, 11, 22, 32, 41, 49, 56, 57, 58, 59, 6, 17, 27,
1914 36, 44, 51, 63, 64, 7, 18, 28, 37, 45, 52, 68, 8,
1915 19, 29, 38, 46, 53, 0, 1, 2, 3, 4, 5, 12, 13,
1916 14, 15, 16, 23, 24, 25, 26, 33, 34, 35, 42, 43, 50}},
1917 {{57, 58, 59, 60, 61, 62, 36, 44, 51, 6, 17, 27, 63,
1918 64, 65, 66, 67, 37, 45, 52, 7, 18, 28, 68, 69, 70,
1919 71, 38, 46, 53, 8, 19, 29, 72, 73, 74, 39, 47, 54,
1920 9, 20, 30, 75, 76, 40, 48, 55, 10, 21, 31, 77, 41,
1921 49, 56, 11, 22, 32, 33, 34, 35, 3, 14, 24, 42, 43,
1922 4, 15, 25, 50, 5, 16, 26, 0, 1, 2, 12, 13, 23},
1923 {72, 73, 74, 60, 65, 69, 39, 47, 54, 9, 20, 30, 75,
1924 76, 61, 66, 70, 40, 48, 55, 10, 21, 31, 77, 62, 67,
1925 71, 41, 49, 56, 11, 22, 32, 57, 58, 59, 36, 44, 51,
1926 6, 17, 27, 63, 64, 37, 45, 52, 7, 18, 28, 68, 38,
1927 46, 53, 8, 19, 29, 33, 34, 35, 3, 14, 24, 42, 43,
1928 4, 15, 25, 50, 5, 16, 26, 0, 1, 2, 12, 13, 23}}}};
1929 s_target = mapDerivIndex2_xxxx[swap_braket][swap_tbra][swap_tket][s];
1933 case BraKet::xs_xx: {
1934 assert(swap_bra ==
false);
1935 assert(swap_braket ==
false);
1936 const unsigned mapDerivIndex2_xsxx[2][45] = {
1937 {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
1938 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
1939 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
1940 36, 37, 38, 39, 40, 41, 42, 43, 44},
1941 {0, 1, 2, 6, 7, 8, 3, 4, 5, 9, 10, 14,
1942 15, 16, 11, 12, 13, 17, 21, 22, 23, 18, 19, 20,
1943 39, 40, 41, 27, 32, 36, 42, 43, 28, 33, 37, 44,
1944 29, 34, 38, 24, 25, 26, 30, 31, 35}};
1945 s_target = mapDerivIndex2_xsxx[swap_tket][s];
1949 case BraKet::xs_xs: {
1950 assert(swap_bra ==
false);
1951 assert(swap_ket ==
false);
1952 assert(swap_braket ==
false);
1958 assert(
false &&
"this backet type not yet supported for 2st geometric derivatives");
1964 "3-rd and higher derivatives not yet generalized");
1967 for (
auto r1 = 0; r1 != nr1; ++r1) {
1968 for (
auto r2 = 0; r2 != nr2; ++r2, src_row_ptr += ncol) {
1969 typedef Eigen::Map<const Matrix> ConstMap;
1970 typedef Eigen::Map<Matrix> Map;
1971 typedef Eigen::Map<Matrix, Eigen::Unaligned,
1972 Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>
1976 ConstMap src_blk_mat(src_row_ptr, nc1, nc2);
1985 const auto tgt_col_idx =
1986 !swap_tket ? r1 * nr2 + r2 : r2 * nr1 + r1;
1987 StridedMap tgt_blk_mat(
1988 tgt_ptr + tgt_col_idx, nr1_tgt, nr2_tgt,
1989 Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>(
1990 nr2_tgt * ncol_tgt, ncol_tgt));
1992 tgt_blk_mat = src_blk_mat.transpose();
1994 tgt_blk_mat = src_blk_mat;
1998 const auto tgt_row_idx =
1999 !swap_tbra ? r1 * nr2 + r2 : r2 * nr1 + r1;
2000 Map tgt_blk_mat(tgt_ptr + tgt_row_idx * ncol, nc1_tgt, nc2_tgt);
2002 tgt_blk_mat = src_blk_mat.transpose();
2004 tgt_blk_mat = src_blk_mat;
2008 std::swap(source, target);
2014 if (source != primdata_[0].targets[s]) {
2015 hotscr += n1234_cart;
2017 assert(set_targets_ &&
"logic error");
2019 targets_[s_target] = source;
2024 assert(set_targets_ &&
"logic error");
2026 targets_[s_target] = source;
2032 for (
auto s = 0; s != ntargets; ++s)
2033 targets_[s] = primdata_[0].targets[s];
2037 #ifdef LIBINT2_ENGINE_TIMERS 2038 const auto t2 = timers.stop(2);
2039 #ifdef LIBINT2_ENGINE_PROFILE_CLASS 2040 class_profiles[id].tform += t2.count();
2045 if (cartesian_shell_normalization() == CartesianShellNormalization::uniform) {
2047 for (
auto s = 0ul; s != targets_.size(); ++s) {
2055 #undef BOOST_PP_NBODY_OPERATOR_LIST 2056 #undef BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE 2057 #undef BOOST_PP_NBODY_OPERATOR_INDEX_LIST 2058 #undef BOOST_PP_NBODY_BRAKET_INDEX_TUPLE 2059 #undef BOOST_PP_NBODY_BRAKET_INDEX_LIST 2060 #undef BOOST_PP_NBODY_DERIV_ORDER_TUPLE 2061 #undef BOOST_PP_NBODY_DERIV_ORDER_LIST 2062 #undef BOOST_PP_NBODYENGINE_MCR3 2063 #undef BOOST_PP_NBODYENGINE_MCR3_ncenter 2064 #undef BOOST_PP_NBODYENGINE_MCR3_default_ncenter 2065 #undef BOOST_PP_NBODYENGINE_MCR3_NCENTER 2066 #undef BOOST_PP_NBODYENGINE_MCR3_OPER 2067 #undef BOOST_PP_NBODYENGINE_MCR3_DERIV 2068 #undef BOOST_PP_NBODYENGINE_MCR3_task 2069 #undef BOOST_PP_NBODYENGINE_MCR3_TASK 2070 #undef BOOST_PP_NBODYENGINE_MCR4 2071 #undef BOOST_PP_NBODYENGINE_MCR5 2072 #undef BOOST_PP_NBODYENGINE_MCR6 2073 #undef BOOST_PP_NBODYENGINE_MCR7 2075 #ifdef LIBINT2_DOES_NOT_INLINE_ENGINE 2076 template any Engine::enforce_params_type<Engine::empty_pod>(
2077 Operator oper,
const Engine::empty_pod& params,
bool throw_if_wrong_type);
2079 template const Engine::target_ptr_vec& Engine::compute<Shell>(
2080 const Shell& first_shell,
const Shell&);
2082 template const Engine::target_ptr_vec& Engine::compute<Shell, Shell>(
2083 const Shell& first_shell,
const Shell&,
const Shell&);
2085 template const Engine::target_ptr_vec& Engine::compute<Shell, Shell, Shell>(
2086 const Shell& first_shell,
const Shell&,
const Shell&,
const Shell&);
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition: shell.h:81
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
svector< Contraction > contr
contractions
Definition: shell.h:105
void uniform_normalize_cartesian_shells(Real *intset, std::array< std::reference_wrapper< const Shell >, N > shells)
rescales cartesian Gaussians to convert from standard to uniform-normalized convention
Definition: cartesian.h:95
bool initialized()
checks if the libint has been initialized.
Definition: initialize.h:61
Array idential to C++0X arrays.
Definition: stdarray_bits.h:14
Iterates over all partitions of a non-negative integer into nonnegative integers in reverse lexicog...
Definition: intpart_iter.h:74
static const Shell & unit()
Definition: shell.h:218
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