19 #ifndef _libint2_src_lib_libint_engine_h_ 20 #define _libint2_src_lib_libint_engine_h_ 22 #include <libint2/cxxstd.h> 23 #if LIBINT2_CPLUSPLUS_STD < 2011 24 # error "libint2/engine.h requires C++11 support" 34 #include <libint2/boys.h> 35 #include <libint2/shell.h> 36 #include <libint2/timer.h> 37 #include <libint2/solidharmonics.h> 38 #include <libint2/any.h> 39 #include <libint2/util/compressed_pair.h> 41 #include <libint2/boost/preprocessor.hpp> 45 #define BOOST_PP_MAKE_TUPLE_INTERNAL(z, i, last) i BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i,last)) 46 #define BOOST_PP_MAKE_TUPLE(n) ( 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 61 #define DEPRECATED __attribute__((deprecated)) 62 #elif defined(_MSC_VER) 63 #define DEPRECATED __declspec(deprecated) 65 #pragma message("WARNING: You need to implement DEPRECATED for this compiler") 71 constexpr
size_t num_geometrical_derivatives(
size_t ncenter,
73 return (deriv_order > 0) ? (num_geometrical_derivatives(ncenter, deriv_order-1) * (3*ncenter+deriv_order-1))/deriv_order : 1;
76 template <
typename T,
unsigned N>
77 typename std::remove_all_extents<T>::type *
78 to_ptr1(T (&a)[N]) {
return reinterpret_cast< 79 typename std::remove_all_extents<T>::type *
82 #if defined(LIBINT2_SUPPORT_ONEBODY) 90 typedef struct {} empty_pod;
107 #define BOOST_PP_ONEBODY_OPERATOR_LIST (overlap, \ 115 #define BOOST_PP_ONEBODY_OPERATOR_INDEX_TUPLE BOOST_PP_MAKE_TUPLE( BOOST_PP_LIST_SIZE(BOOST_PP_ONEBODY_OPERATOR_LIST) ) 116 #define BOOST_PP_ONEBODY_OPERATOR_INDEX_LIST BOOST_PP_TUPLE_TO_LIST( BOOST_PP_ONEBODY_OPERATOR_INDEX_TUPLE ) 119 #define BOOST_PP_ONEBODY_DERIV_ORDER_TUPLE BOOST_PP_MAKE_TUPLE( BOOST_PP_INC(LIBINT2_DERIV_ONEBODY_ORDER) ) 120 #define BOOST_PP_ONEBODY_DERIV_ORDER_LIST BOOST_PP_TUPLE_TO_LIST( BOOST_PP_ONEBODY_DERIV_ORDER_TUPLE ) 134 OneBodyEngine() : type_(_invalid), primdata_(), stack_size_(0), lmax_(-1) {}
145 template <
typename Params = empty_pod>
150 Params params = empty_pod()) :
152 primdata_(max_nprim * max_nprim),
154 deriv_order_(deriv_order),
155 params_(enforce_params_type(type,params)),
156 fm_eval_(type ==
nuclear ? coulomb_core_eval_t::instance(2*max_l+deriv_order, 1e-25) : decltype(fm_eval_){})
165 primdata_(
std::move(other.primdata_)),
166 stack_size_(other.stack_size_),
168 hard_lmax_(other.hard_lmax_),
169 deriv_order_(other.deriv_order_),
170 params_(
std::move(other.params_)),
171 fm_eval_(
std::move(other.fm_eval_)),
172 scratch_(
std::move(other.scratch_)),
173 scratch2_(other.scratch2_),
174 buildfnptrs_(other.buildfnptrs_) {
180 primdata_(other.primdata_.size()),
181 stack_size_(other.stack_size_),
183 deriv_order_(other.deriv_order_),
184 params_(other.params_),
185 fm_eval_(other.fm_eval_) {
196 primdata_ = std::move(other.primdata_);
197 stack_size_ = other.stack_size_;
199 hard_lmax_ = other.hard_lmax_;
200 deriv_order_ = other.deriv_order_;
201 params_ = std::move(other.params_);
202 fm_eval_ = std::move(other.fm_eval_);
203 scratch_ = std::move(other.scratch_);
204 scratch2_ = other.scratch2_;
205 buildfnptrs_ = other.buildfnptrs_;
212 primdata_.resize(other.primdata_.size());
213 stack_size_ = other.stack_size_;
215 deriv_order_ = other.deriv_order_;
216 params_ = other.params_;
217 fm_eval_ = other.fm_eval_;
224 template <
typename Params>
231 template <
typename Params>
232 DEPRECATED
void set_q(
const Params& params) {
241 const unsigned int num_operator_geometrical_derivatives = (type_ ==
nuclear) ? this->nparams() : 0;
242 const auto ncenters = 2 + num_operator_geometrical_derivatives;
243 return nopers() * num_geometrical_derivatives(ncenters,deriv_order_);
255 unsigned int deriv_order,
256 unsigned int ncenters,
257 unsigned int omitted_center) {
258 auto ncenters_reduced = ncenters - 1;
259 auto nderiv_1d = 3 * ncenters_reduced;
260 assert(deriv_idx < num_geometrical_derivatives(ncenters_reduced,deriv_order));
261 switch (deriv_order) {
262 case 1:
return deriv_idx >= omitted_center*3 ? deriv_idx+3 : deriv_idx;
263 default: assert(deriv_order<=1);
275 assert(s1.ncontr() == 1 && s2.ncontr() == 1);
277 const auto l1 = s1.
contr[0].l;
278 const auto l2 = s2.
contr[0].l;
281 if (type_ ==
nuclear and nparams() == 0)
282 throw std::runtime_error(
"libint2::OneBodyEngine<nuclear>, but no charges found; forgot to call set_params()?");
284 const auto n1 = s1.size();
285 const auto n2 = s2.size();
286 const auto n12 = n1 * n2;
287 const auto ncart1 = s1.cartesian_size();
288 const auto ncart2 = s2.cartesian_size();
289 const auto ncart12 = ncart1 * ncart2;
290 const auto nops = nopers();
292 const auto tform_to_solids = s1.
contr[0].pure || s2.
contr[0].pure;
295 const auto nprim1 = s1.nprim();
296 const auto nprim2 = s2.nprim();
297 const auto nprimpairs = nprim1 * nprim2;
298 assert(nprimpairs <= primdata_.size());
300 auto nparam_sets = nparams();
305 const auto geometry_independent_operator = type_ ==
overlap || type_ ==
kinetic;
306 const auto num_deriv_centers_computed = geometry_independent_operator ? 1 : 2;
307 auto num_shellsets_computed = nopers() *
308 num_geometrical_derivatives(num_deriv_centers_computed,
311 const auto target_buf_size = num_shellsets_computed * ncart12;
317 const auto accumulate_ints_in_scratch = (type_ ==
nuclear);
320 const auto lmax = std::max(l1, l2);
321 assert (lmax <= lmax_);
325 auto cartesian_ints = primdata_[0].stack;
328 const auto compute_directly = lmax == 0 && deriv_order_ == 0 && (type_ ==
overlap || type_ ==
nuclear);
329 if (compute_directly) {
330 primdata_[0].stack[0] = 0;
332 else if (accumulate_ints_in_scratch)
333 memset(static_cast<void*>(&scratch_[0]), 0,
sizeof(real_t)*target_buf_size);
336 for(
auto pset=0u; pset!=nparam_sets; ++pset) {
338 if (type_!=
nuclear) assert(nparam_sets == 1);
341 for(
auto p1=0; p1!=nprim1; ++p1) {
342 for(
auto p2=0; p2!=nprim2; ++p2, ++p12) {
343 compute_primdata(primdata_[p12],s1,s2,p1,p2,pset);
346 primdata_[0].contrdepth = p12;
348 if (compute_directly) {
349 auto& result = cartesian_ints[0];
352 for(
auto p12=0; p12 != primdata_[0].contrdepth; ++p12)
353 result += primdata_[p12]._0_Overlap_0_x[0]
354 * primdata_[p12]._0_Overlap_0_y[0]
355 * primdata_[p12]._0_Overlap_0_z[0];
358 for(
auto p12=0; p12 != primdata_[0].contrdepth; ++p12)
359 result += primdata_[p12].LIBINT_T_S_ELECPOT_S(0)[0];
364 primdata_[0].targets[0] = cartesian_ints;
368 buildfnptrs_[s1.
contr[0].l*hard_lmax_ + s2.
contr[0].l](&primdata_[0]);
369 cartesian_ints = primdata_[0].targets[0];
371 if (accumulate_ints_in_scratch) {
372 cartesian_ints = &scratch_[0];
373 std::transform(primdata_[0].targets[0], primdata_[0].targets[0] + target_buf_size,
375 &scratch_[0], std::plus<real_t>());
378 if (deriv_order_ > 0){
379 const auto nints_per_center = target_buf_size/2;
382 auto dest = &scratch_[0] + (2+pset)*nints_per_center;
383 auto src = primdata_[0].targets[0];
384 for(
auto i=0; i!=nints_per_center; ++i) {
387 src = primdata_[0].targets[0] + nints_per_center;
388 for(
auto i=0; i!=nints_per_center; ++i) {
391 num_shellsets_computed+=3;
399 auto result = cartesian_ints;
400 if (tform_to_solids) {
402 auto* spherical_ints = (cartesian_ints == &scratch_[0]) ? scratch2_ : &scratch_[0];
403 result = spherical_ints;
407 for(
auto s=0ul; s!=num_shellsets_computed; ++s, cartesian_ints+=ncart12) {
418 spherical_ints = result + n12 * s_target;
420 libint2::solidharmonics::tform(l1, l2, cartesian_ints, spherical_ints);
423 if (s1.
contr[0].pure)
424 libint2::solidharmonics::tform_rows(l1, n2, cartesian_ints, spherical_ints);
426 libint2::solidharmonics::tform_cols(n1, l2, cartesian_ints, spherical_ints);
434 if (deriv_order_ > 0 && geometry_independent_operator) {
435 assert(deriv_order_ == 1);
437 const auto nints_computed = n12*num_shellsets_computed;
441 if (not tform_to_solids) {
442 const auto stack_size_remaining = stack_size_ - (result-primdata_[0].stack) - nints_computed;
443 const auto copy_to_scratch2 = stack_size_remaining < nints_computed;
444 if (copy_to_scratch2) {
447 std::copy(result, result + nints_computed, scratch2_);
452 const auto src = result;
453 const auto dest = result + nints_computed;
454 for(
auto f=0ul; f!=nints_computed; ++f) {
462 inline void compute_primdata(
Libint_t& primdata,
464 size_t p1,
size_t p2,
469 std::vector<Libint_t> primdata_;
475 std::shared_ptr<coulomb_core_eval_t> fm_eval_;
476 std::vector<real_t> scratch_;
479 typedef void (*buildfnptr_t)(
const Libint_t*);
480 buildfnptr_t* buildfnptrs_;
482 void reset_scratch() {
483 const auto ncart_max = (lmax_+1)*(lmax_+2)/2;
484 const auto target_shellset_size =
nshellsets() * ncart_max * ncart_max;
489 const auto need_extra_large_scratch = stack_size_ < target_shellset_size;
490 scratch_.resize(need_extra_large_scratch ? 2*target_shellset_size : target_shellset_size);
491 scratch2_ = need_extra_large_scratch ? &scratch_[target_shellset_size] : primdata_[0].stack;
495 assert(libint2::initialized());
496 assert(deriv_order_ <= LIBINT2_DERIV_ONEBODY_ORDER);
498 #define BOOST_PP_ONEBODYENGINE_MCR2(r,product) \ 499 if (type_ == BOOST_PP_TUPLE_ELEM(2,0,product) && deriv_order_ == BOOST_PP_TUPLE_ELEM(2,1,product) ) {\ 500 assert(lmax_ <= BOOST_PP_CAT(LIBINT2_MAX_AM_ , \ 501 BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \ 502 BOOST_PP_TUPLE_ELEM(2,0,product) ) \ 505 LIBINT2_PREFIXED_NAME( BOOST_PP_CAT( \ 506 BOOST_PP_CAT(libint2_need_memory_ , \ 507 BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \ 508 BOOST_PP_TUPLE_ELEM(2,0,product) ) \ 510 BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \ 511 BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \ 515 LIBINT2_PREFIXED_NAME( BOOST_PP_CAT( \ 516 BOOST_PP_CAT(libint2_init_ , \ 517 BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \ 518 BOOST_PP_TUPLE_ELEM(2,0,product) ) \ 520 BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \ 521 BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \ 524 )(&primdata_[0], lmax_, 0); \ 525 buildfnptrs_ = to_ptr1( \ 526 LIBINT2_PREFIXED_NAME( BOOST_PP_CAT( \ 527 BOOST_PP_CAT(libint2_build_ , \ 528 BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \ 529 BOOST_PP_TUPLE_ELEM(2,0,product) ) \ 531 BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \ 532 BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \ 537 hard_lmax_ = BOOST_PP_CAT( \ 540 BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \ 541 BOOST_PP_TUPLE_ELEM(2,0,product) ), \ 542 BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \ 543 BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \ 551 BOOST_PP_LIST_FOR_EACH_PRODUCT ( BOOST_PP_ONEBODYENGINE_MCR2, 2, (BOOST_PP_ONEBODY_OPERATOR_INDEX_LIST, BOOST_PP_ONEBODY_DERIV_ORDER_LIST) )
557 if (primdata_.size() != 0) {
559 #define BOOST_PP_ONEBODYENGINE_MCR3(r,product) \ 560 LIBINT2_PREFIXED_NAME( BOOST_PP_CAT( \ 561 BOOST_PP_CAT(libint2_cleanup_ , \ 562 BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \ 563 BOOST_PP_TUPLE_ELEM(2,0,product) ) \ 565 BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \ 566 BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \ 572 BOOST_PP_LIST_FOR_EACH_PRODUCT ( BOOST_PP_ONEBODYENGINE_MCR3, 2, (BOOST_PP_ONEBODY_OPERATOR_INDEX_LIST, BOOST_PP_ONEBODY_DERIV_ORDER_LIST) )
580 inline unsigned int nparams()
const;
581 inline unsigned int nopers()
const;
585 template <
typename Params>
587 const Params& params,
588 bool throw_if_wrong_type = not std::is_same<Params,empty_pod>::value);
594 static oper_params_type default_params() {
return oper_params_type{}; }
595 static constexpr
unsigned int nopers = 1;
601 static oper_params_type default_params() {
return oper_params_type{}; }
602 static constexpr
unsigned int nopers = 1;
608 static oper_params_type default_params() {
return oper_params_type{{0.0,0.0,0.0}}; }
609 static constexpr
unsigned int nopers = 4;
624 inline unsigned int OneBodyEngine::nparams()
const {
633 inline unsigned int OneBodyEngine::nopers()
const {
635 #define BOOST_PP_ONEBODYENGINE_MCR4(r,data,i,elem) case i : return operator_traits< static_cast<OneBodyEngine::operator_type> ( i ) >::nopers; 636 BOOST_PP_LIST_FOR_EACH_I ( BOOST_PP_ONEBODYENGINE_MCR4, _, BOOST_PP_ONEBODY_OPERATOR_LIST)
642 template <
typename Params>
644 const Params& params,
645 bool throw_if_wrong_type) {
649 #define BOOST_PP_ONEBODYENGINE_MCR5(r,data,i,elem) \ 651 if (std::is_same<Params,operator_traits< static_cast<operator_type> ( i ) >::oper_params_type>::value) \ 654 if (throw_if_wrong_type) throw std::bad_cast(); \ 655 result = operator_traits<static_cast<operator_type> ( i ) >::default_params(); \ 659 BOOST_PP_LIST_FOR_EACH_I ( BOOST_PP_ONEBODYENGINE_MCR5, _, BOOST_PP_ONEBODY_OPERATOR_LIST)
667 inline void OneBodyEngine::compute_primdata(
Libint_t& primdata,
669 size_t p1,
size_t p2,
672 const auto& A = s1.
O;
673 const auto& B = s2.
O;
675 const auto alpha1 = s1.
alpha[p1];
676 const auto alpha2 = s2.
alpha[p2];
678 const auto c1 = s1.
contr[0].coeff[p1];
679 const auto c2 = s2.
contr[0].coeff[p2];
681 const auto gammap = alpha1 + alpha2;
682 const auto oogammap = 1.0 / gammap;
683 const auto rhop_over_alpha1 = alpha2 * oogammap;
684 const auto rhop = alpha1 * rhop_over_alpha1;
685 const auto Px = (alpha1 * A[0] + alpha2 * B[0]) * oogammap;
686 const auto Py = (alpha1 * A[1] + alpha2 * B[1]) * oogammap;
687 const auto Pz = (alpha1 * A[2] + alpha2 * B[2]) * oogammap;
688 const auto AB_x = A[0] - B[0];
689 const auto AB_y = A[1] - B[1];
690 const auto AB_z = A[2] - B[2];
691 const auto AB2_x = AB_x * AB_x;
692 const auto AB2_y = AB_y * AB_y;
693 const auto AB2_z = AB_z * AB_z;
694 const auto AB2 = AB2_x + AB2_y + AB2_z;
696 assert (LIBINT2_SHELLQUARTET_SET == LIBINT2_SHELLQUARTET_SET_STANDARD);
700 #if LIBINT2_DEFINED(eri,PA_x) 701 primdata.PA_x[0] = Px - A[0];
703 #if LIBINT2_DEFINED(eri,PA_y) 704 primdata.PA_y[0] = Py - A[1];
706 #if LIBINT2_DEFINED(eri,PA_z) 707 primdata.PA_z[0] = Pz - A[2];
712 #if LIBINT2_DEFINED(eri,PB_x) 713 primdata.PB_x[0] = Px - B[0];
715 #if LIBINT2_DEFINED(eri,PB_y) 716 primdata.PB_y[0] = Py - B[1];
718 #if LIBINT2_DEFINED(eri,PB_z) 719 primdata.PB_z[0] = Pz - B[2];
725 #if LIBINT2_DEFINED(eri,BO_x) 726 primdata.BO_x[0] = B[0] - O[0];
728 #if LIBINT2_DEFINED(eri,BO_y) 729 primdata.BO_y[0] = B[1] - O[1];
731 #if LIBINT2_DEFINED(eri,BO_z) 732 primdata.BO_z[0] = B[2] - O[2];
736 #if LIBINT2_DEFINED(eri,oo2z) 737 primdata.oo2z[0] = 0.5*oogammap;
742 const auto& C = params[oset].second;
743 #if LIBINT2_DEFINED(eri,PC_x) 744 primdata.PC_x[0] = Px - C[0];
746 #if LIBINT2_DEFINED(eri,PC_y) 747 primdata.PC_y[0] = Py - C[1];
749 #if LIBINT2_DEFINED(eri,PC_z) 750 primdata.PC_z[0] = Pz - C[2];
753 #if LIBINT2_DEFINED(eri,AB_x) 754 primdata.AB_x[0] = A[0] - B[0];
756 #if LIBINT2_DEFINED(eri,AB_y) 757 primdata.AB_y[0] = A[1] - B[1];
759 #if LIBINT2_DEFINED(eri,AB_z) 760 primdata.AB_z[0] = A[2] - B[2];
765 decltype(c1) sqrt_PI(1.77245385090551602729816748334);
766 const auto xyz_pfac = sqrt_PI * sqrt(oogammap);
767 const auto ovlp_ss_x = exp(- rhop * AB2_x) * xyz_pfac * c1 * c2;
768 const auto ovlp_ss_y = exp(- rhop * AB2_y) * xyz_pfac;
769 const auto ovlp_ss_z = exp(- rhop * AB2_z) * xyz_pfac;
771 primdata._0_Overlap_0_x[0] = ovlp_ss_x;
772 primdata._0_Overlap_0_y[0] = ovlp_ss_y;
773 primdata._0_Overlap_0_z[0] = ovlp_ss_z;
775 if (type_ ==
kinetic || (deriv_order_ > 0)) {
776 #if LIBINT2_DEFINED(eri,two_alpha0_bra) 777 primdata.two_alpha0_bra[0] = 2.0 * alpha1;
779 #if LIBINT2_DEFINED(eri,two_alpha0_ket) 780 primdata.two_alpha0_ket[0] = 2.0 * alpha2;
785 #if LIBINT2_DEFINED(eri,rho12_over_alpha1) || LIBINT2_DEFINED(eri,rho12_over_alpha2) 786 if (deriv_order_ > 0) {
787 #if LIBINT2_DEFINED(eri,rho12_over_alpha1) 788 primdata.rho12_over_alpha1[0] = rhop_over_alpha1;
790 #if LIBINT2_DEFINED(eri,rho12_over_alpha2) 791 primdata.rho12_over_alpha2[0] = alpha1 * oogammap;
795 #if LIBINT2_DEFINED(eri,PC_x) && LIBINT2_DEFINED(eri,PC_y) && LIBINT2_DEFINED(eri,PC_z) 796 const auto PC2 = primdata.PC_x[0] * primdata.PC_x[0] +
797 primdata.PC_y[0] * primdata.PC_y[0] +
798 primdata.PC_z[0] * primdata.PC_z[0];
799 const auto U = gammap * PC2;
800 const auto mmax = s1.
contr[0].l + s2.
contr[0].l + deriv_order_;
801 auto* fm_ptr = &(primdata.LIBINT_T_S_ELECPOT_S(0)[0]);
802 fm_eval_->eval(fm_ptr, U, mmax);
804 decltype(U) two_o_sqrt_PI(1.12837916709551257389615890312);
806 const auto pfac = - q * sqrt(gammap) * two_o_sqrt_PI * ovlp_ss_x * ovlp_ss_y * ovlp_ss_z;
807 const auto m_fence = mmax + 1;
808 for(
auto m=0; m!=m_fence; ++m) {
816 #undef BOOST_PP_ONEBODY_OPERATOR_LIST 817 #undef BOOST_PP_ONEBODY_OPERATOR_INDEX_TUPLE 818 #undef BOOST_PP_ONEBODY_OPERATOR_INDEX_LIST 819 #undef BOOST_PP_ONEBODY_DERIV_ORDER_TUPLE 820 #undef BOOST_PP_ONEBODY_DERIV_ORDER_LIST 821 #undef BOOST_PP_ONEBODYENGINE_MCR2 822 #undef BOOST_PP_ONEBODYENGINE_MCR3 823 #undef BOOST_PP_ONEBODYENGINE_MCR4 824 #undef BOOST_PP_ONEBODYENGINE_MCR5 826 #endif // LIBINT2_SUPPORT_ONEBODY 862 typedef struct {} oper_params_type;
866 typedef ContractedGaussianGeminal oper_params_type;
870 typedef ContractedGaussianGeminal oper_params_type;
874 typedef ContractedGaussianGeminal oper_params_type;
878 void operator()(
double* Gm,
882 const static auto one_over_two_pi = 1.0 / (2.0 * M_PI);
883 const auto G0 = exp(-T) * rho * one_over_two_pi;
884 std::fill(Gm, Gm+mmax+1, G0);
888 typedef libint2::GenericGmEval<delta_gm_eval> core_eval_type;
889 typedef struct {} oper_params_type;
892 #ifdef LIBINT2_SUPPORT_ERI 899 template <MultiplicativeSphericalTwoBodyKernel Kernel>
906 set_precision(std::numeric_limits<real_t>::epsilon());
928 real_t precision = std::numeric_limits<real_t>::epsilon(),
929 const oper_params_type& oparams = oper_params_type()) :
930 primdata_(max_nprim * max_nprim * max_nprim * max_nprim),
931 spbra_(max_nprim), spket_(max_nprim),
932 lmax_(max_l), deriv_order_(deriv_order),
933 core_eval_pack_(core_eval_type::instance(4*lmax_ + deriv_order,
934 std::numeric_limits<real_t>::epsilon()),
938 set_precision(precision);
940 init_core_ints_params(oparams);
946 primdata_(
std::move(other.primdata_)),
947 spbra_(
std::move(other.spbra_)), spket_(
std::move(other.spket_)),
948 lmax_(other.lmax_), deriv_order_(other.deriv_order_),
949 precision_(other.precision_), ln_precision_(other.ln_precision_),
950 core_eval_pack_(
std::move(other.core_eval_pack_)),
951 core_ints_params_(
std::move(other.core_ints_params_)),
952 scratch_(
std::move(other.scratch_)) {
957 primdata_(other.primdata_),
958 spbra_(other.spbra_), spket_(other.spket_),
959 lmax_(other.lmax_), deriv_order_(other.deriv_order_),
960 precision_(other.precision_), ln_precision_(other.ln_precision_),
961 core_eval_pack_(other.core_eval_pack_), core_ints_params_(other.core_ints_params_)
973 primdata_ = std::move(other.primdata_);
975 deriv_order_ = other.deriv_order_;
976 precision_ = other.precision_;
977 ln_precision_ = other.ln_precision_;
978 core_eval_pack_ = std::move(other.core_eval_pack_);
979 core_ints_params_ = std::move(other.core_ints_params_);
980 scratch_ = std::move(other.scratch_);
987 primdata_ = other.primdata_;
989 deriv_order_ = other.deriv_order_;
990 precision_ = other.precision_;
991 ln_precision_ = other.ln_precision_;
992 core_eval_pack_ = other.core_eval_pack_;
993 core_ints_params_ = other.core_ints_params_;
998 static bool skip_core_ints;
1000 #ifdef LIBINT2_ENGINE_TIMERS 1004 # ifdef LIBINT2_ENGINE_PROFILE_CLASS 1007 template <
typename I>
class_id(I l0, I l1, I l2, I l3) {
1013 bool operator<(
const class_id& other)
const {
1014 return ordinal(l) < ordinal(other.l);
1016 static size_t ordinal(
const size_t (&l)[4]) {
1017 return ((l[0]*LIBINT2_MAX_AM_ERI + l[1])*LIBINT2_MAX_AM_ERI + l[2])*LIBINT2_MAX_AM_ERI + l[3];
1020 std::ostringstream oss;
1021 oss <<
"(" << Shell::am_symbol(l[0]) << Shell::am_symbol(l[1])
1022 <<
"|" << Shell::am_symbol(l[2]) << Shell::am_symbol(l[3]) <<
")";
1036 prereqs = build_hrr = build_vrr = tform = 0.;
1037 nprimset = nshellset = 0;
1040 std::map<class_id, class_profile> class_profiles;
1056 assert(tbra1.ncontr() == 1 && tbra2.ncontr() == 1 &&
1057 tket1.ncontr() == 1 && tket2.ncontr() == 1);
1060 assert(deriv_order_ == 0);
1062 #if LIBINT2_SHELLQUARTET_SET == LIBINT2_SHELLQUARTET_SET_STANDARD // standard angular momentum ordering 1063 auto swap_bra = (tbra1.
contr[0].l < tbra2.
contr[0].l);
1064 auto swap_ket = (tket1.
contr[0].l < tket2.
contr[0].l);
1065 auto swap_braket = (tbra1.
contr[0].l + tbra2.
contr[0].l > tket1.
contr[0].l + tket2.
contr[0].l);
1066 #else // orca angular momentum ordering 1067 auto swap_bra = (tbra1.
contr[0].l > tbra2.
contr[0].l);
1068 auto swap_ket = (tket1.
contr[0].l > tket2.
contr[0].l);
1069 auto swap_braket = (tbra1.
contr[0].l + tbra2.
contr[0].l < tket1.
contr[0].l + tket2.
contr[0].l);
1071 const auto& bra1 = swap_braket ? (swap_ket ? tket2 : tket1) : (swap_bra ? tbra2 : tbra1);
1072 const auto& bra2 = swap_braket ? (swap_ket ? tket1 : tket2) : (swap_bra ? tbra1 : tbra2);
1073 const auto& ket1 = swap_braket ? (swap_bra ? tbra2 : tbra1) : (swap_ket ? tket2 : tket1);
1074 const auto& ket2 = swap_braket ? (swap_bra ? tbra1 : tbra2) : (swap_ket ? tket1 : tket2);
1076 const bool tform = tbra1.contr[0].pure || tbra2.contr[0].pure || tket1.contr[0].pure || tket2.contr[0].pure;
1077 const bool use_scratch = (swap_braket || swap_bra || swap_ket || tform);
1080 auto nprim_bra1 = bra1.nprim();
1081 auto nprim_bra2 = bra2.nprim();
1082 auto nprim_ket1 = ket1.nprim();
1083 auto nprim_ket2 = ket2.nprim();
1086 auto lmax = std::max(std::max(bra1.contr[0].l, bra2.contr[0].l), std::max(ket1.contr[0].l, ket2.contr[0].l));
1087 assert (lmax <= lmax_);
1089 primdata_[0].stack[0] = 0;
1091 #ifdef LIBINT2_ENGINE_PROFILE_CLASS 1092 class_id id(bra1.contr[0].l, bra2.contr[0].l, ket1.contr[0].l, ket2.contr[0].l);
1093 if (class_profiles.find(
id) == class_profiles.end()) {
1095 class_profiles[id] = dummy;
1100 #ifdef LIBINT2_ENGINE_TIMERS 1108 spbra_.init(bra1, bra2, ln_precision_);
1109 spket_.init(ket1, ket2, ln_precision_);
1110 const auto npbra = spbra_.primpairs.size();
1111 const auto npket = spket_.primpairs.size();
1112 for(
auto pb=0; pb!=npbra; ++pb) {
1113 for(
auto pk=0; pk!=npket; ++pk) {
1114 if (spbra_.primpairs[pb].scr + spket_.primpairs[pk].scr > ln_precision_) {
1115 if (compute_primdata(primdata_[p],bra1,bra2,ket1,ket2,spbra_,pb,spket_,pk)) {
1121 primdata_[0].contrdepth = p;
1124 #ifdef LIBINT2_ENGINE_TIMERS 1125 const auto t0 = timers.
stop(0);
1126 # ifdef LIBINT2_ENGINE_PROFILE_CLASS 1127 class_profiles[id].prereqs += t0.count();
1128 if (primdata_[0].contrdepth != 0) {
1129 class_profiles[id].nshellset += 1;
1130 class_profiles[id].nprimset += primdata_[0].contrdepth;
1136 if (primdata_[0].contrdepth == 0) {
1137 const size_t n = bra1.size() * bra2.size() * ket1.size() * ket2.size();
1138 memset(primdata_[0].stack, 0,
sizeof(real_t)*n);
1139 return primdata_[0].stack;
1142 real_t* result =
nullptr;
1145 #ifdef LIBINT2_ENGINE_TIMERS 1148 auto& stack = primdata_[0].stack[0];
1150 for(
auto p=0; p != primdata_[0].contrdepth; ++p)
1151 stack += primdata_[p].LIBINT_T_SS_EREP_SS(0)[0];
1152 primdata_[0].targets[0] = primdata_[0].stack;
1153 #ifdef LIBINT2_ENGINE_TIMERS 1154 const auto t1 = timers.
stop(1);
1155 # ifdef LIBINT2_ENGINE_PROFILE_CLASS 1156 class_profiles[id].build_vrr += t1.count();
1160 result = primdata_[0].targets[0];
1163 #ifdef LIBINT2_ENGINE_TIMERS 1164 # ifdef LIBINT2_PROFILE 1165 const auto t1_hrr_start = primdata_[0].timers->read(0);
1166 const auto t1_vrr_start = primdata_[0].timers->read(1);
1170 LIBINT2_PREFIXED_NAME(libint2_build_eri)[bra1.contr[0].l][bra2.contr[0].l][ket1.contr[0].l][ket2.contr[0].l](&primdata_[0]);
1171 #ifdef LIBINT2_ENGINE_TIMERS 1172 const auto t1 = timers.
stop(1);
1173 # ifdef LIBINT2_ENGINE_PROFILE_CLASS 1174 # ifndef LIBINT2_PROFILE 1175 class_profiles[id].build_vrr += t1.count();
1177 class_profiles[id].build_hrr += primdata_[0].timers->read(0) - t1_hrr_start;
1178 class_profiles[id].build_vrr += primdata_[0].timers->read(1) - t1_vrr_start;
1183 result = primdata_[0].targets[0];
1185 #ifdef LIBINT2_ENGINE_TIMERS 1192 constexpr
auto using_scalar_real = std::is_same<double,real_t>::value || std::is_same<float,real_t>::value;
1193 static_assert(using_scalar_real,
"Libint2 C++11 API only supports fundamental real types");
1194 typedef Eigen::Matrix<real_t, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > Matrix;
1197 const auto nr1_cart = bra1.cartesian_size();
1198 const auto nr2_cart = bra2.cartesian_size();
1199 const auto nc1_cart = ket1.cartesian_size();
1200 const auto nc2_cart = ket2.cartesian_size();
1201 const auto ncol_cart = nc1_cart * nc2_cart;
1202 const auto nr1 = bra1.size();
1203 const auto nr2 = bra2.size();
1204 const auto nc1 = ket1.size();
1205 const auto nc2 = ket2.size();
1206 const auto nrow = nr1 * nr2;
1207 const auto ncol = nc1 * nc2;
1210 const auto nr1_tgt = tbra1.size();
1211 const auto nr2_tgt = tbra2.size();
1212 const auto nc1_tgt = tket1.size();
1213 const auto nc2_tgt = tket2.size();
1214 const auto ncol_tgt = nc1_tgt * nc2_tgt;
1217 auto mainbuf = result;
1218 auto scratchbuf = &scratch_[0];
1219 if (bra1.contr[0].pure) {
1220 libint2::solidharmonics::transform_first(bra1.contr[0].l, nr2_cart*ncol_cart,
1221 mainbuf, scratchbuf);
1222 std::swap(mainbuf, scratchbuf);
1224 if (bra2.contr[0].pure) {
1225 libint2::solidharmonics::transform_inner(bra1.size(), bra2.contr[0].l, ncol_cart,
1226 mainbuf, scratchbuf);
1227 std::swap(mainbuf, scratchbuf);
1229 if (ket1.contr[0].pure) {
1230 libint2::solidharmonics::transform_inner(nrow, ket1.contr[0].l, nc2_cart,
1231 mainbuf, scratchbuf);
1232 std::swap(mainbuf, scratchbuf);
1234 if (ket2.contr[0].pure) {
1235 libint2::solidharmonics::transform_last(bra1.size()*bra2.size()*ket1.size(), ket2.contr[0].l,
1236 mainbuf, scratchbuf);
1237 std::swap(mainbuf, scratchbuf);
1241 const auto* src_row_ptr = mainbuf;
1242 auto tgt_ptr = scratchbuf;
1243 for(
auto r1=0; r1!=nr1; ++r1) {
1244 for(
auto r2=0; r2!=nr2; ++r2, src_row_ptr+=ncol) {
1246 typedef Eigen::Map<const Matrix> ConstMap;
1247 typedef Eigen::Map<Matrix> Map;
1248 typedef Eigen::Map<Matrix, Eigen::Unaligned, Eigen::Stride<Eigen::Dynamic,Eigen::Dynamic> > StridedMap;
1251 ConstMap src_blk_mat(src_row_ptr, nc1, nc2);
1257 const auto tgt_col_idx = !swap_ket ? r1 * nr2 + r2 : r2 * nr1 + r1;
1258 StridedMap tgt_blk_mat(tgt_ptr + tgt_col_idx, nr1_tgt, nr2_tgt,
1259 Eigen::Stride<Eigen::Dynamic,Eigen::Dynamic>(nr2_tgt*ncol_tgt,ncol_tgt));
1261 tgt_blk_mat = src_blk_mat.transpose();
1263 tgt_blk_mat = src_blk_mat;
1267 const auto tgt_row_idx = !swap_bra ? r1 * nr2 + r2 : r2 * nr1 + r1;
1268 Map tgt_blk_mat(tgt_ptr + tgt_row_idx*ncol, nc1_tgt, nc2_tgt);
1270 tgt_blk_mat = src_blk_mat.transpose();
1272 tgt_blk_mat = src_blk_mat;
1278 result = scratchbuf;
1282 #ifdef LIBINT2_ENGINE_TIMERS 1283 const auto t2 = timers.
stop(2);
1284 # ifdef LIBINT2_ENGINE_PROFILE_CLASS 1285 class_profiles[id].tform += t2.count();
1305 ln_precision_ = std::numeric_limits<real_t>::lowest();
1309 ln_precision_ = std::log(precision_);
1318 void print_timers() {
1319 #ifdef LIBINT2_ENGINE_TIMERS 1320 std::cout <<
"timers: prereq = " << timers.
read(0);
1321 # ifndef LIBINT2_PROFILE // if libint's profiling was on, engine's build timer will include its overhead 1323 std::cout <<
" build = " << timers.
read(1);
1325 std::cout <<
" tform = " << timers.
read(2) << std::endl;
1327 #ifdef LIBINT2_PROFILE 1328 std::cout <<
"build timers: hrr = " << primdata_[0].timers->read(0)
1329 <<
" vrr = " << primdata_[0].timers->read(1) << std::endl;
1331 #ifdef LIBINT2_ENGINE_TIMERS 1332 # ifdef LIBINT2_ENGINE_PROFILE_CLASS 1333 for(
const auto& p: class_profiles) {
1334 printf(
"{\"%s\", %10.5lf, %10.5lf, %10.5lf, %10.5lf, %ld, %ld},\n",
1335 p.first.to_string().c_str(),
1349 inline bool compute_primdata(
Libint_t& primdata,
1357 std::vector<Libint_t> primdata_;
1360 size_t deriv_order_;
1362 real_t ln_precision_;
1367 core_eval_scratch_type>
1369 core_eval_pack_type core_eval_pack_;
1371 typedef oper_params_type core_ints_params_type;
1372 core_ints_params_type core_ints_params_;
1374 void init_core_ints_params(
const oper_params_type& oparams);
1376 std::vector<real_t> scratch_;
1381 assert(libint2::initialized());
1382 assert(lmax_ <= LIBINT2_MAX_AM_ERI);
1383 assert(deriv_order_ <= LIBINT2_DERIV_ERI_ORDER);
1385 const auto ncart_max = (lmax_+1)*(lmax_+2)/2;
1386 const auto max_shellpair_size = ncart_max * ncart_max;
1387 const auto max_shellset_size = max_shellpair_size * max_shellpair_size;
1389 switch(deriv_order_) {
1392 libint2_init_eri(&primdata_[0], lmax_, 0);
1393 scratch_.resize(max_shellset_size);
1396 #if LIBINT2_DERIV_ERI_ORDER > 0 1397 libint2_init_eri1(&primdata_[0], lmax_, 0);
1398 scratch_.resize(9 * max_shellset_size);
1402 #if LIBINT2_DERIV_ERI_ORDER > 1 1403 libint2_init_eri2(&primdata_[0], lmax_, 0);
1404 scratch_.resize(45 * max_shellset_size);
1407 default: assert(deriv_order_ < 3);
1410 #ifdef LIBINT2_ENGINE_TIMERS 1413 #ifdef LIBINT2_PROFILE 1414 primdata_[0].timers->set_now_overhead(25);
1419 if (primdata_.size() != 0) {
1420 switch(deriv_order_) {
1423 libint2_cleanup_eri(&primdata_[0]);
1426 #if LIBINT2_DERIV_ERI_ORDER > 0 1427 libint2_cleanup_eri1(&primdata_[0]);
1431 #if LIBINT2_DERIV_ERI_ORDER > 1 1432 libint2_cleanup_eri2(&primdata_[0]);
1448 engine->core_eval_pack_.first()->eval(Fm, T, mmax);
1459 engine->core_eval_pack_.first()->eval(Gm, rho, T, mmax, engine->core_ints_params_, &engine->core_eval_pack_.second());
1469 engine->core_eval_pack_.first()->eval(Gm, rho, T, mmax, engine->core_ints_params_);
1479 engine->core_eval_pack_.first()->eval(Gm, rho, T, mmax, engine->core_ints_params_);
1489 engine->core_eval_pack_.first()->eval(Gm, rho, T, mmax);
1494 template <MultiplicativeSphericalTwoBodyKernel Kernel>
1497 template <MultiplicativeSphericalTwoBodyKernel Kernel>
1506 const auto& A = sbra1.
O;
1507 const auto& B = sbra2.
O;
1508 const auto& C = sket1.
O;
1509 const auto& D = sket2.
O;
1510 const auto& AB = spbra.AB;
1511 const auto& CD = spket.AB;
1513 const auto& spbrapp = spbra.primpairs[pbra];
1514 const auto& spketpp = spket.primpairs[pket];
1515 const auto& pbra1 = spbrapp.p1;
1516 const auto& pbra2 = spbrapp.p2;
1517 const auto& pket1 = spketpp.p1;
1518 const auto& pket2 = spketpp.p2;
1520 const auto alpha0 = sbra1.
alpha[pbra1];
1521 const auto alpha1 = sbra2.
alpha[pbra2];
1522 const auto alpha2 = sket1.
alpha[pket1];
1523 const auto alpha3 = sket2.
alpha[pket2];
1525 const auto c0 = sbra1.
contr[0].coeff[pbra1];
1526 const auto c1 = sbra2.
contr[0].coeff[pbra2];
1527 const auto c2 = sket1.
contr[0].coeff[pket1];
1528 const auto c3 = sket2.
contr[0].coeff[pket2];
1530 const auto amtot = sbra1.
contr[0].l + sket1.
contr[0].l +
1533 const auto gammap = alpha0 + alpha1;
1534 const auto oogammap = spbrapp.one_over_gamma;
1535 const auto rhop = alpha0 * alpha1 * oogammap;
1537 const auto gammaq = alpha2 + alpha3;
1538 const auto oogammaq = spketpp.one_over_gamma;
1539 const auto rhoq = alpha2 * alpha3 * oogammaq;
1541 const auto& P = spbrapp.P;
1542 const auto& Q = spketpp.P;
1543 const auto PQx = P[0] - Q[0];
1544 const auto PQy = P[1] - Q[1];
1545 const auto PQz = P[2] - Q[2];
1546 const auto PQ2 = PQx * PQx + PQy * PQy + PQz * PQz;
1548 const auto K12 = spbrapp.K * spketpp.K;
1549 decltype(K12) two_times_M_PI_to_25(34.986836655249725693);
1550 const auto gammapq = gammap + gammaq;
1551 const auto sqrt_gammapq = sqrt(gammapq);
1552 const auto oogammapq = 1.0 / (gammapq);
1553 auto pfac = two_times_M_PI_to_25 * K12 * sqrt_gammapq * oogammapq;
1554 pfac *= c0 * c1 * c2 * c3;
1556 if (std::abs(pfac) < precision_)
1559 const auto rho = gammap * gammaq * oogammapq;
1560 const auto T = PQ2*rho;
1561 auto* fm_ptr = &(primdata.LIBINT_T_SS_EREP_SS(0)[0]);
1562 const auto mmax = amtot + deriv_order_;
1564 if (!skip_core_ints)
1567 for(
auto m=0; m!=mmax+1; ++m) {
1574 #if LIBINT2_DEFINED(eri,PA_x) 1575 primdata.PA_x[0] = P[0] - A[0];
1577 #if LIBINT2_DEFINED(eri,PA_y) 1578 primdata.PA_y[0] = P[1] - A[1];
1580 #if LIBINT2_DEFINED(eri,PA_z) 1581 primdata.PA_z[0] = P[2] - A[2];
1583 #if LIBINT2_DEFINED(eri,PB_x) 1584 primdata.PB_x[0] = P[0] - B[0];
1586 #if LIBINT2_DEFINED(eri,PB_y) 1587 primdata.PB_y[0] = P[1] - B[1];
1589 #if LIBINT2_DEFINED(eri,PB_z) 1590 primdata.PB_z[0] = P[2] - B[2];
1593 #if LIBINT2_DEFINED(eri,QC_x) 1594 primdata.QC_x[0] = Q[0] - C[0];
1596 #if LIBINT2_DEFINED(eri,QC_y) 1597 primdata.QC_y[0] = Q[1] - C[1];
1599 #if LIBINT2_DEFINED(eri,QC_z) 1600 primdata.QC_z[0] = Q[2] - C[2];
1602 #if LIBINT2_DEFINED(eri,QD_x) 1603 primdata.QD_x[0] = Q[0] - D[0];
1605 #if LIBINT2_DEFINED(eri,QD_y) 1606 primdata.QD_y[0] = Q[1] - D[1];
1608 #if LIBINT2_DEFINED(eri,QD_z) 1609 primdata.QD_z[0] = Q[2] - D[2];
1612 #if LIBINT2_DEFINED(eri,AB_x) 1613 primdata.AB_x[0] = AB[0];
1615 #if LIBINT2_DEFINED(eri,AB_y) 1616 primdata.AB_y[0] = AB[1];
1618 #if LIBINT2_DEFINED(eri,AB_z) 1619 primdata.AB_z[0] = AB[2];
1621 #if LIBINT2_DEFINED(eri,BA_x) 1622 primdata.BA_x[0] = -AB[0];
1624 #if LIBINT2_DEFINED(eri,BA_y) 1625 primdata.BA_y[0] = -AB[1];
1627 #if LIBINT2_DEFINED(eri,BA_z) 1628 primdata.BA_z[0] = -AB[2];
1631 #if LIBINT2_DEFINED(eri,CD_x) 1632 primdata.CD_x[0] = CD[0];
1634 #if LIBINT2_DEFINED(eri,CD_y) 1635 primdata.CD_y[0] = CD[1];
1637 #if LIBINT2_DEFINED(eri,CD_z) 1638 primdata.CD_z[0] = CD[2];
1640 #if LIBINT2_DEFINED(eri,DC_x) 1641 primdata.DC_x[0] = -CD[0];
1643 #if LIBINT2_DEFINED(eri,DC_y) 1644 primdata.DC_y[0] = -CD[1];
1646 #if LIBINT2_DEFINED(eri,DC_z) 1647 primdata.DC_z[0] = -CD[2];
1650 const auto gammap_o_gammapgammaq = oogammapq * gammap;
1651 const auto gammaq_o_gammapgammaq = oogammapq * gammaq;
1653 const auto Wx = (gammap_o_gammapgammaq * P[0] + gammaq_o_gammapgammaq * Q[0]);
1654 const auto Wy = (gammap_o_gammapgammaq * P[1] + gammaq_o_gammapgammaq * Q[1]);
1655 const auto Wz = (gammap_o_gammapgammaq * P[2] + gammaq_o_gammapgammaq * Q[2]);
1657 #if LIBINT2_DEFINED(eri,WP_x) 1658 primdata.WP_x[0] = Wx - P[0];
1660 #if LIBINT2_DEFINED(eri,WP_y) 1661 primdata.WP_y[0] = Wy - P[1];
1663 #if LIBINT2_DEFINED(eri,WP_z) 1664 primdata.WP_z[0] = Wz - P[2];
1666 #if LIBINT2_DEFINED(eri,WQ_x) 1667 primdata.WQ_x[0] = Wx - Q[0];
1669 #if LIBINT2_DEFINED(eri,WQ_y) 1670 primdata.WQ_y[0] = Wy - Q[1];
1672 #if LIBINT2_DEFINED(eri,WQ_z) 1673 primdata.WQ_z[0] = Wz - Q[2];
1675 #if LIBINT2_DEFINED(eri,oo2z) 1676 primdata.oo2z[0] = 0.5*oogammap;
1678 #if LIBINT2_DEFINED(eri,oo2e) 1679 primdata.oo2e[0] = 0.5*oogammaq;
1681 #if LIBINT2_DEFINED(eri,oo2ze) 1682 primdata.oo2ze[0] = 0.5*oogammapq;
1684 #if LIBINT2_DEFINED(eri,roz) 1685 primdata.roz[0] = rho*oogammap;
1687 #if LIBINT2_DEFINED(eri,roe) 1688 primdata.roe[0] = rho*oogammaq;
1692 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x) 1693 primdata.TwoPRepITR_pfac0_0_0_x[0] = - (alpha1*AB[0] + alpha3*CD[0]) * oogammap;
1695 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y) 1696 primdata.TwoPRepITR_pfac0_0_0_y[0] = - (alpha1*AB[1] + alpha3*CD[1]) * oogammap;
1698 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z) 1699 primdata.TwoPRepITR_pfac0_0_0_z[0] = - (alpha1*AB[2] + alpha3*CD[2]) * oogammap;
1701 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x) 1702 primdata.TwoPRepITR_pfac0_1_0_x[0] = - (alpha1*AB[0] + alpha3*CD[0]) * oogammaq;
1704 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y) 1705 primdata.TwoPRepITR_pfac0_1_0_y[0] = - (alpha1*AB[1] + alpha3*CD[1]) * oogammaq;
1707 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z) 1708 primdata.TwoPRepITR_pfac0_1_0_z[0] = - (alpha1*AB[2] + alpha3*CD[2]) * oogammaq;
1710 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x) 1711 primdata.TwoPRepITR_pfac0_0_1_x[0] = (alpha0*AB[0] + alpha2*CD[0]) * oogammap;
1713 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y) 1714 primdata.TwoPRepITR_pfac0_0_1_y[0] = (alpha0*AB[1] + alpha2*CD[1]) * oogammap;
1716 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z) 1717 primdata.TwoPRepITR_pfac0_0_1_z[0] = (alpha0*AB[2] + alpha2*CD[2]) * oogammap;
1719 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x) 1720 primdata.TwoPRepITR_pfac0_1_1_x[0] = (alpha0*AB[0] + alpha2*CD[0]) * oogammaq;
1722 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y) 1723 primdata.TwoPRepITR_pfac0_1_1_y[0] = (alpha0*AB[1] + alpha2*CD[1]) * oogammaq;
1725 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z) 1726 primdata.TwoPRepITR_pfac0_1_1_z[0] = (alpha0*AB[2] + alpha2*CD[2]) * oogammaq;
1728 #if LIBINT2_DEFINED(eri,eoz) 1729 primdata.eoz[0] = gammaq*oogammap;
1731 #if LIBINT2_DEFINED(eri,zoe) 1732 primdata.zoe[0] = gammap*oogammaq;
1736 if (deriv_order_ > 0) {
1737 #if LIBINT2_DEFINED(eri,alpha1_rho_over_zeta2) 1738 primdata.alpha1_rho_over_zeta2[0] = alpha0 * (oogammap * gammaq_o_gammapgammaq);
1740 #if LIBINT2_DEFINED(eri,alpha2_rho_over_zeta2) 1741 primdata.alpha2_rho_over_zeta2[0] = alpha1 * (oogammap * gammaq_o_gammapgammaq);
1743 #if LIBINT2_DEFINED(eri,alpha3_rho_over_eta2) 1744 primdata.alpha3_rho_over_eta2[0] = alpha2 * (oogammaq * gammap_o_gammapgammaq);
1746 #if LIBINT2_DEFINED(eri,alpha4_rho_over_eta2) 1747 primdata.alpha4_rho_over_eta2[0] = alpha3 * (oogammaq * gammap_o_gammapgammaq);
1749 #if LIBINT2_DEFINED(eri,alpha1_over_zetapluseta) 1750 primdata.alpha1_over_zetapluseta[0] = alpha0 * oogammapq;
1752 #if LIBINT2_DEFINED(eri,alpha2_over_zetapluseta) 1753 primdata.alpha2_over_zetapluseta[0] = alpha1 * oogammapq;
1755 #if LIBINT2_DEFINED(eri,alpha3_over_zetapluseta) 1756 primdata.alpha3_over_zetapluseta[0] = alpha2 * oogammapq;
1758 #if LIBINT2_DEFINED(eri,alpha4_over_zetapluseta) 1759 primdata.alpha4_over_zetapluseta[0] = alpha3 * oogammapq;
1761 #if LIBINT2_DEFINED(eri,rho12_over_alpha1) 1762 primdata.rho12_over_alpha1[0] = alpha1 * oogammap;
1764 #if LIBINT2_DEFINED(eri,rho12_over_alpha2) 1765 primdata.rho12_over_alpha2[0] = alpha0 * oogammap;
1767 #if LIBINT2_DEFINED(eri,rho34_over_alpha3) 1768 primdata.rho34_over_alpha3[0] = alpha3 * oogammaq;
1770 #if LIBINT2_DEFINED(eri,rho34_over_alpha4) 1771 primdata.rho34_over_alpha4[0] = alpha2 * oogammaq;
1773 #if LIBINT2_DEFINED(eri,two_alpha0_bra) 1774 primdata.two_alpha0_bra[0] = 2.0 * alpha0;
1776 #if LIBINT2_DEFINED(eri,two_alpha0_ket) 1777 primdata.two_alpha0_ket[0] = 2.0 * alpha1;
1779 #if LIBINT2_DEFINED(eri,two_alpha1_bra) 1780 primdata.two_alpha1_bra[0] = 2.0 * alpha2;
1782 #if LIBINT2_DEFINED(eri,two_alpha1_ket) 1783 primdata.two_alpha1_ket[0] = 2.0 * alpha3;
1792 const oper_params_type& oparams) {
1796 const auto ng = oparams.size();
1797 core_ints_params_.reserve(ng*(ng+1)/2);
1798 for(
size_t b=0; b<ng; ++b)
1799 for(
size_t k=0; k<=b; ++k) {
1800 const auto gexp = oparams[b].first + oparams[k].first;
1801 const auto gcoeff = oparams[b].second * oparams[k].second * (b == k ? 1 : 2);
1802 const auto gcoeff_rescaled = 4 * oparams[b].first * oparams[k].first * gcoeff;
1803 core_ints_params_.push_back(std::make_pair(gexp, gcoeff_rescaled));
1807 template <MultiplicativeSphericalTwoBodyKernel Kernel>
1809 const oper_params_type& oparams) {
1810 core_ints_params_ = oparams;
1814 #endif // LIBINT2_SUPPORT_ERI std::vector< real_t > alpha
exponents
Definition: shell.h:100
ShellPair pre-computes shell-pair data, primitive pairs are screened to finite precision.
Definition: shell.h:302
DEPRECATED typedef operator_type integral_type
alias to operator_type for backward compatibility to pre-05/13/2015 code
Definition: engine.h:124
TwoBodyEngine & operator=(TwoBodyEngine &&other)
move assignment
Definition: engine.h:972
void set_params(const Params ¶ms)
resets operator parameters; this may be useful if need to compute Coulomb potential integrals over ba...
Definition: engine.h:225
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition: shell.h:78
TwoBodyEngine(size_t max_nprim, int max_l, int deriv_order=0, real_t precision=std::numeric_limits< real_t >::epsilon(), const oper_params_type &oparams=oper_params_type())
Constructs a (usable) TwoBodyEngine.
Definition: engine.h:926
std::array< double, 3 > oper_params_type
Cartesian coordinates of the origin with respect to which the dipole moment is defined.
Definition: engine.h:607
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
( cGTG) ( cGTG)
Definition: engine.h:833
Definition: stdarray.h:18
TwoBodyEngine & operator=(const TwoBodyEngine &other)
(deep) copy assignment
Definition: engine.h:985
unsigned int nshellsets() const
reports the number of shell sets that each call to compute() produces.
Definition: engine.h:240
some evaluators need thread-local scratch, but most don't
Definition: boys.h:1325
std::vector< std::pair< double, double > > ContractedGaussianGeminal
contracted Gaussian geminal = , represented as a vector of { , } pairs
Definition: engine.h:839
emulates boost::any
Definition: any.h:17
const real_t * compute(const libint2::Shell &tbra1, const libint2::Shell &tbra2, const libint2::Shell &tket1, const libint2::Shell &tket2)
computes shell set of integrals
Definition: engine.h:1046
describes operator sets given by OneBodyOperator
Definition: engine.h:128
OneBodyEngine::operator_traits< OneBodyEngine::emultipole1 >::oper_params_type oper_params_type
Cartesian coordinates of the origin with respect to which the multipole moments are defined...
Definition: engine.h:619
Coulomb potential due to point charges.
Definition: engine.h:99
overlap + (Cartesian) electric dipole moment, , where is relative to origin
Definition: engine.h:100
static unsigned int to_target_deriv_index(unsigned int deriv_idx, unsigned int deriv_order, unsigned int ncenters, unsigned int omitted_center)
Given the Cartesian geometric derivative index that refers to center set (0...n-1) with one center om...
Definition: engine.h:254
OneBodyEngine & operator=(const OneBodyEngine &other)
(deep) copy assignment
Definition: engine.h:210
std::array< real_t, 3 > O
origin
Definition: shell.h:102
contracted Gaussian geminal times Coulomb
Definition: engine.h:832
std::string to_string(const T &x)
Converts x to its string representation.
Definition: entity.h:71
Definition: compressed_pair.h:27
Definition: engine.h:1005
void set_precision(real_t prec)
this specifies target precision for computing the integrals.
Definition: engine.h:1302
dur_t stop(size_t t)
stops timer t
Definition: timer.h:63
OneBodyEngine & operator=(OneBodyEngine &&other)
move assignment is default
Definition: engine.h:194
TwoBodyEngine(TwoBodyEngine &&other)
move constructor
Definition: engine.h:945
overlap
Definition: engine.h:97
emultipole2 + (Cartesian) electric octupole moment,
Definition: engine.h:102
const real_t * compute(const libint2::Shell &s1, const libint2::Shell &s2)
computes shell set of integrals
Definition: engine.h:271
TwoBodyEngine(const TwoBodyEngine &other)
(deep) copy constructor
Definition: engine.h:956
double read(size_t t) const
reads value (in seconds) of timer t , converted to double
Definition: timer.h:70
Definition: engine.h:1026
Definition: test_eri_rys.cc:46
Computes the Boys function, $ F_m (T) = ^1 u^{2m} (-T u^2) \, { d}u $, using 3-rd order Chebyshev int...
Definition: boys.h:239
Array idential to C++0X arrays.
Definition: stdarray_bits.h:14
electronic kinetic energy, i.e.
Definition: engine.h:98
OneBodyEngine(operator_type type, size_t max_nprim, int max_l, int deriv_order=0, Params params=empty_pod())
Constructs a (usable) OneBodyEngine.
Definition: engine.h:146
emultipole1 + (Cartesian) electric quadrupole moment,
Definition: engine.h:101
MultiplicativeSphericalTwoBodyKernel
types of multiplicative spherically-symmetric two-body kernels known by TwoBodyEngine ...
Definition: engine.h:829
void start(size_t t)
starts timer t
Definition: timer.h:58
operator_type
types of operators (operator sets) supported by OneBodyEngine.
Definition: engine.h:96
contracted Gaussian geminal =
Definition: engine.h:831
DEPRECATED void set_q(const Params ¶ms)
alias to set_params() for backward compatibility with pre-05/13/2015 code
Definition: engine.h:232
OneBodyEngine::operator_traits< OneBodyEngine::emultipole1 >::oper_params_type oper_params_type
Cartesian coordinates of the origin with respect to which the multipole moments are defined...
Definition: engine.h:613
TwoBodyEngine()
creates a default (unusable) TwoBodyEngine
Definition: engine.h:905
OneBodyEngine(OneBodyEngine &&other)
move constructor
Definition: engine.h:163
std::vector< std::pair< double, std::array< double, 3 > > > oper_params_type
point charges and their positions
Definition: engine.h:600
OneBodyEngine()
creates a default (unusable) OneBodyEngine; to be used as placeholder for copying a usable engine ...
Definition: engine.h:134
void set_now_overhead(size_t ns)
use this to report the overhead of now() call; if set, the reported timings will be adjusted for this...
Definition: timer.h:53
TwoBodyEngine computes (ab|O|cd) (i.e.
Definition: engine.h:900
Computes the Boys function, $ F_m (T) = ^1 u^{2m} (-T u^2) \, { d}u $, using Taylor interpolation of ...
Definition: boys.h:752
OneBodyEngine(const OneBodyEngine &other)
(deep) copy constructor
Definition: engine.h:178
real_t precision() const
Definition: engine.h:1314
std::vector< Contraction > contr
contractions
Definition: shell.h:101
OneBodyEngine computes integrals of operators (or operator sets) given by OneBodyOperator::operator_t...
Definition: engine.h:87