20 #ifndef _libint2_src_bin_libint_vrr11r12kg1211_h_ 21 #define _libint2_src_bin_libint_vrr11r12kg1211_h_ 23 #include <generic_rr.h> 24 #include <r12kg12_11_11.h> 34 template <
class BFSet,
int part, FunctionPosition where>
37 GenIntegralSet_11_11<BFSet,R12kG12,mType> >
45 static const unsigned int max_nchildren = 8;
47 using ParentType::Instance;
50 static bool directional() {
return ParentType::default_directional(); }
53 using ParentType::RecurrenceRelation::expr_;
54 using ParentType::RecurrenceRelation::nflops_;
55 using ParentType::target_;
56 using ParentType::is_simple;
61 static std::string descr() {
return "VRR"; }
65 std::string generate_label()
const 67 typedef typename TargetType::AuxIndexType
mType;
68 static SafePtr<mType> aux0(
new mType(0u));
70 os << descr() <<
"P" << part <<
to_string(where)
71 << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
75 #if LIBINT_ENABLE_GENERIC_CODE 76 bool has_generic(
const SafePtr<CompilationParameters>& cparams)
const;
79 std::string generic_header()
const;
81 std::string generic_instance(
const SafePtr<CodeContext>& context,
const SafePtr<CodeSymbols>& args)
const;
85 template <
class F,
int part, FunctionPosition where>
92 const int K = Tint->oper()->descr().K();
94 const R12kG12 oK((R12kG12Descr(K)));
95 const unsigned int m = Tint->aux()->elem(0);
96 const F _1 = unit<F>(dir);
104 if (a.contracted() ||
108 Tint->oper()->descr().contracted())
117 const bool deriv = dA.
zero() ==
false ||
118 dB.
zero() ==
false ||
119 dC.
zero() ==
false ||
121 assert(deriv ==
false);
131 if (part == 0 && where == InBra) {
132 F a(Tint->bra(0,0) - _1);
138 auto ABCD_m = factory.
make_child(a,b,c,d,m,oK);
139 auto ABCD_mp1 = factory.
make_child(a,b,c,d,m+1,oK);
140 if (is_simple()) { expr_ = Vector(
"PA")[dir] * ABCD_m + Vector(
"WP")[dir] * ABCD_mp1; nflops_+=3; }
142 const F& am1 = a - _1;
144 auto Am1BCD_m = factory.
make_child(am1,b,c,d,m,oK);
145 auto Am1BCD_mp1 = factory.
make_child(am1,b,c,d,m+1,oK);
146 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2z") * (Am1BCD_m - Scalar(
"roz") * Am1BCD_mp1); nflops_+=5; }
148 const F& bm1 = b - _1;
150 auto ABm1CD_m = factory.
make_child(a,bm1,c,d,m,oK);
151 auto ABm1CD_mp1 = factory.
make_child(a,bm1,c,d,m+1,oK);
152 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2z") * (ABm1CD_m - Scalar(
"roz") * ABm1CD_mp1); nflops_+=5; }
154 const F& cm1 = c - _1;
156 auto ABCm1D_mp1 = factory.
make_child(a,b,cm1,d,m+1,oK);
157 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2ze") * ABCm1D_mp1; nflops_+=3; }
159 const F& dm1 = d - _1;
161 auto ABCDm1_mp1 = factory.
make_child(a,b,c,dm1,m+1,oK);
162 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2ze") * ABCDm1_mp1; nflops_+=3; }
167 if (part == 0 && where == InKet) {
169 F b(Tint->ket(0,0) - _1);
174 auto ABCD_m = factory.
make_child(a,b,c,d,m,oK);
175 auto ABCD_mp1 = factory.
make_child(a,b,c,d,m+1,oK);
176 if (is_simple()) { expr_ = Vector(
"PB")[dir] * ABCD_m + Vector(
"WP")[dir] * ABCD_mp1; nflops_+=3; }
178 const F& am1 = a - _1;
180 auto Am1BCD_m = factory.
make_child(am1,b,c,d,m,oK);
181 auto Am1BCD_mp1 = factory.
make_child(am1,b,c,d,m+1,oK);
182 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2z") * (Am1BCD_m - Scalar(
"roz") * Am1BCD_mp1); nflops_+=5; }
184 const F& bm1 = b - _1;
186 auto ABm1CD_m = factory.
make_child(a,bm1,c,d,m,oK);
187 auto ABm1CD_mp1 = factory.
make_child(a,bm1,c,d,m+1,oK);
188 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2z") * (ABm1CD_m - Scalar(
"roz") * ABm1CD_mp1); nflops_+=5; }
190 const F& cm1 = c - _1;
192 auto ABCm1D_mp1 = factory.
make_child(a,b,cm1,d,m+1,oK);
193 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2ze") * ABCm1D_mp1; nflops_+=3; }
195 const F& dm1 = d - _1;
197 auto ABCDm1_mp1 = factory.
make_child(a,b,c,dm1,m+1,oK);
198 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2ze") * ABCDm1_mp1; nflops_+=3; }
203 if (part == 1 && where == InBra) {
206 F c(Tint->bra(1,0) - _1);
210 auto ABCD_m = factory.
make_child(a,b,c,d,m,oK);
211 auto ABCD_mp1 = factory.
make_child(a,b,c,d,m+1,oK);
212 if (is_simple()) { expr_ = Vector(
"QC")[dir] * ABCD_m + Vector(
"WQ")[dir] * ABCD_mp1; nflops_+=3; }
214 const F& cm1 = c - _1;
216 auto ABCm1D_m = factory.
make_child(a,b,cm1,d,m,oK);
217 auto ABCm1D_mp1 = factory.
make_child(a,b,cm1,d,m+1,oK);
218 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2e") * (ABCm1D_m - Scalar(
"roe") * ABCm1D_mp1); nflops_+=5; }
220 const F& dm1 = d - _1;
222 auto ABCDm1_m = factory.
make_child(a,b,c,dm1,m,oK);
223 auto ABCDm1_mp1 = factory.
make_child(a,b,c,dm1,m+1,oK);
224 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2e") * (ABCDm1_m - Scalar(
"roe") * ABCDm1_mp1); nflops_+=5; }
226 const F& am1 = a - _1;
228 auto Am1BCD_mp1 = factory.
make_child(am1,b,c,d,m+1,oK);
229 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2ze") * Am1BCD_mp1; nflops_+=3; }
231 const F& bm1 = b - _1;
233 auto ABm1CD_mp1 = factory.
make_child(a,bm1,c,d,m+1,oK);
234 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2ze") * ABm1CD_mp1; nflops_+=3; }
239 if (part == 1 && where == InKet) {
243 F d(Tint->ket(1,0) - _1);
246 auto ABCD_m = factory.
make_child(a,b,c,d,m,oK);
247 auto ABCD_mp1 = factory.
make_child(a,b,c,d,m+1,oK);
248 if (is_simple()) { expr_ = Vector(
"QD")[dir] * ABCD_m + Vector(
"WQ")[dir] * ABCD_mp1; nflops_+=3; }
250 const F& cm1 = c - _1;
252 auto ABCm1D_m = factory.
make_child(a,b,cm1,d,m,oK);
253 auto ABCm1D_mp1 = factory.
make_child(a,b,cm1,d,m+1,oK);
254 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2e") * (ABCm1D_m - Scalar(
"roe") * ABCm1D_mp1); nflops_+=5; }
256 const F& dm1 = d - _1;
258 auto ABCDm1_m = factory.
make_child(a,b,c,dm1,m,oK);
259 auto ABCDm1_mp1 = factory.
make_child(a,b,c,dm1,m+1,oK);
260 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2e") * (ABCDm1_m - Scalar(
"roe") * ABCDm1_mp1); nflops_+=5; }
262 const F& am1 = a - _1;
264 auto Am1BCD_mp1 = factory.
make_child(am1,b,c,d,m+1,oK);
265 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2ze") * Am1BCD_mp1; nflops_+=3; }
267 const F& bm1 = b - _1;
269 auto ABm1CD_mp1 = factory.
make_child(a,bm1,c,d,m+1,oK);
270 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2ze") * ABm1CD_mp1; nflops_+=3; }
279 throw std::logic_error(
"VRR_11_R12kG12_11<I,F,K,part,where>::children_and_expr_Kge0() -- nonzero auxiliary quantum detected.");
285 F b(Tint->ket(0,0) - _1);
286 F d(Tint->ket(1,0) - _1);
291 F a(Tint->bra(0,0) - _1);
292 F c(Tint->bra(1,0) - _1);
305 if (part == 0 && where == InBra) {
306 F a(Tint->bra(0,0) - _1);
312 auto ABCD_K = factory.
make_child(a,b,c,d,0u,oK);
313 if (is_simple()) { expr_ = Vector(
"R12kG12_pfac0_0")[dir] * ABCD_K; nflops_+=1; }
314 const F& am1 = a - _1;
316 auto Am1BCD_K = factory.
make_child(am1,b,c,d,0u,oK);
317 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"R12kG12_pfac1_0") * Am1BCD_K; nflops_+=3; }
319 const F& cm1 = c - _1;
321 auto ABCm1D_K = factory.
make_child(a,b,cm1,d,0u,oK);
322 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"R12kG12_pfac2") * ABCm1D_K; nflops_+=3; }
325 const R12kG12 oKm2((R12kG12Descr(K-2)));
326 auto Ap1BCD_Km2 = factory.
make_child(a+_1,b,c,d,0u,oKm2);
327 auto ABCp1D_Km2 = factory.
make_child(a,b,c+_1,d,0u,oKm2);
328 auto ABCD_Km2 = factory.
make_child(a,b,c,d,0u,oKm2);
329 if (is_simple()) { expr_ += Scalar((
double)K) * Scalar(
"R12kG12_pfac3_0")
330 * (Ap1BCD_Km2 - ABCp1D_Km2 + Vector(
"R12kG12_pfac4_0")[dir] * ABCD_Km2); nflops_+=6; }
335 if (part == 1 && where == InBra) {
338 F c(Tint->bra(1,0) - _1);
342 auto ABCD_K = factory.
make_child(a,b,c,d,0u,oK);
343 if (is_simple()) { expr_ = Vector(
"R12kG12_pfac0_1")[dir] * ABCD_K; nflops_+=1; }
344 const F& cm1 = c - _1;
346 auto ABCm1D_K = factory.
make_child(a,b,cm1,d,0u,oK);
347 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"R12kG12_pfac1_1") * ABCm1D_K; nflops_+=3; }
349 const F& am1 = a - _1;
351 auto Am1BCD_K = factory.
make_child(am1,b,c,d,0u,oK);
352 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"R12kG12_pfac2") * Am1BCD_K; nflops_+=3; }
355 const R12kG12 oKm2((R12kG12Descr(K-2)));
356 auto ABCp1D_Km2 = factory.
make_child(a,b,c+_1,d,0u,oKm2);
357 auto Ap1BCD_Km2 = factory.
make_child(a+_1,b,c,d,0u,oKm2);
358 auto ABCD_Km2 = factory.
make_child(a,b,c,d,0u,oKm2);
359 if (is_simple()) { expr_ += Scalar((
double)K) * Scalar(
"R12kG12_pfac3_1")
360 * (ABCp1D_Km2 - Ap1BCD_Km2 + Vector(
"R12kG12_pfac4_1")[dir] * ABCD_Km2); nflops_+=6; }
368 if (part == 0 && where == InKet) {
370 F b(Tint->ket(0,0) - _1);
375 auto ABCD_K = factory.
make_child(a,b,c,d,0u,oK);
376 if (is_simple()) { expr_ = Vector(
"R12kG12_pfac0_0")[dir] * ABCD_K; nflops_+=1; }
377 const F& bm1 = b - _1;
379 auto ABm1CD_K = factory.
make_child(a,bm1,c,d,0u,oK);
380 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"R12kG12_pfac1_0") * ABm1CD_K; nflops_+=3; }
382 const F& dm1 = d - _1;
384 auto ABCDm1_K = factory.
make_child(a,b,c,dm1,0u,oK);
385 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"R12kG12_pfac2") * ABCDm1_K; nflops_+=3; }
388 const R12kG12 oKm2((R12kG12Descr(K-2)));
389 auto ABp1CD_Km2 = factory.
make_child(a,b+_1,c,d,0u,oKm2);
390 auto ABCDp1_Km2 = factory.
make_child(a,b,c,d+_1,0u,oKm2);
391 auto ABCD_Km2 = factory.
make_child(a,b,c,d,0u,oKm2);
392 if (is_simple()) { expr_ += Scalar((
double)K) * Scalar(
"R12kG12_pfac3_0")
393 * (ABp1CD_Km2 - ABCDp1_Km2 + Vector(
"R12kG12_pfac4_0")[dir] * ABCD_Km2); nflops_+=6; }
398 if (part == 1 && where == InKet) {
402 F d(Tint->ket(1,0) - _1);
405 auto ABCD_K = factory.
make_child(a,b,c,d,0u,oK);
406 if (is_simple()) { expr_ = Vector(
"R12kG12_pfac0_1")[dir] * ABCD_K; nflops_+=1; }
407 const F& dm1 = d - _1;
409 auto ABCDm1_K = factory.
make_child(a,b,c,dm1,0u,oK);
410 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"R12kG12_pfac1_1") * ABCDm1_K; nflops_+=3; }
412 const F& bm1 = b - _1;
414 auto ABm1CD_K = factory.
make_child(a,bm1,c,d,0u,oK);
415 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"R12kG12_pfac2") * ABm1CD_K; nflops_+=3; }
418 const R12kG12 oKm2((R12kG12Descr(K-2)));
419 auto ABCDp1_Km2 = factory.
make_child(a,b,c,d+_1,0u,oKm2);
420 auto ABp1CD_Km2 = factory.
make_child(a,b+_1,c,d,0u,oKm2);
421 auto ABCD_Km2 = factory.
make_child(a,b,c,d,0u,oKm2);
422 if (is_simple()) { expr_ += Scalar((
double)K) * Scalar(
"R12kG12_pfac3_1")
423 * (ABCDp1_Km2 - ABp1CD_Km2 + Vector(
"R12kG12_pfac4_1")[dir] * ABCD_Km2); nflops_+=6; }
433 #if LIBINT_ENABLE_GENERIC_CODE 434 template <
class F,
int part, FunctionPosition where>
438 F sh_a(target_->bra(0,0));
439 F sh_b(target_->ket(0,0));
440 F sh_c(target_->bra(1,0));
441 F sh_d(target_->ket(1,0));
442 const unsigned int max_opt_am = cparams->max_am_opt();
446 sh_b.zero() && sh_d.zero() &&
447 (sh_a.norm() > std::max(2*max_opt_am,1u) ||
448 sh_c.norm() > std::max(2*max_opt_am,1u)
450 (sh_a.norm() > 1u && sh_c.norm() > 1u)
454 sh_a.zero() && sh_c.zero() &&
455 (sh_b.norm() > std::max(2*max_opt_am,1u) ||
456 sh_d.norm() > std::max(2*max_opt_am,1u)
458 (sh_b.norm() > 1u && sh_d.norm() > 1u)
464 template <
class F,
int part, FunctionPosition where>
468 F sh_a(target_->bra(0,0));
469 F sh_b(target_->ket(0,0));
470 F sh_c(target_->bra(1,0));
471 F sh_d(target_->ket(1,0));
472 const bool xsxs = sh_b.zero() && sh_d.zero();
473 const bool sxsx = sh_a.zero() && sh_c.zero();
474 const int K = target_->oper()->descr().K();
477 return std::string(
"OSVRR_xs_xs.h");
479 return std::string(
"OSVRR_sx_sx.h");
482 return std::string(
"VRR_r12kg12_xs_xs.h");
487 template <
class F,
int part, FunctionPosition where>
491 const int K = target_->oper()->descr().K();
492 std::ostringstream oss;
493 F sh_a(target_->bra(0,0));
494 F sh_b(target_->ket(0,0));
495 F sh_c(target_->bra(1,0));
496 F sh_d(target_->ket(1,0));
497 const bool xsxs = sh_b.zero() && sh_d.zero();
498 const bool sxsx = sh_a.zero() && sh_c.zero();
500 bool ahlrichs_simplification =
false;
503 unit_s = sh_b.is_unit();
506 unit_s = sh_a.is_unit();
509 oss <<
"using namespace libint2;" << endl;
512 oss <<
"libint2::OSVRR_xs_xs<" << part <<
"," << sh_a.norm() <<
"," << sh_c.norm() <<
",";
513 if (not ahlrichs_simplification) oss << (unit_s ?
"true," :
"false,");
514 oss << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
517 oss <<
"libint2::OSVRR_sx_sx<" << part <<
"," << sh_b.norm() <<
"," << sh_d.norm() <<
",";
518 if (not ahlrichs_simplification) oss << (unit_s ?
"true," :
"false,");
519 oss << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
524 oss <<
"libint2::VRR_r12kg12_xs_xs<" << part <<
"," << sh_a.norm() <<
"," << sh_c.norm() <<
",";
525 oss << K <<
"," << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
529 oss <<
"libint2::VRR_r12kg12_xs_xs<" << part <<
"," << sh_b.norm() <<
"," << sh_d.norm() <<
",";
530 oss << K <<
"," << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
533 oss <<
">::compute(inteval";
537 oss <<
"," << args->symbol(0);
538 if (not ahlrichs_simplification && unit_s)
540 const unsigned int nargs = args->n();
541 for(
unsigned int a=1; a<nargs; a++) {
542 oss <<
"," << args->symbol(a);
546 const unsigned int nargs = args->n();
547 for(
unsigned int a=0; a<nargs; a++) {
548 oss <<
"," << args->symbol(a);
553 for(
unsigned int a=0; a<3; ++a) {
554 oss <<
",(LIBINT2_REALTYPE*)0";
562 #endif // #if !LIBINT_ENABLE_GENERIC_CODE R12_k_G12 is a two-body operator of form r_{12}^k * exp(- * r_{12}), where k is an integer and is a ...
Definition: oper.h:298
const SafePtr< DGVertex > & make_child(const F &A, const F &B, const F &C, const F &D, const AuxIndexType &aux=AuxIndexType(), const OperType &oper=OperType())
make_child
Definition: generic_rr.h:157
TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition: bfset.h:879
static bool directional()
Default directionality.
Definition: vrr_11_r12kg12_11.h:50
VRR Recurrence Relation for 2-e integrals of the R12_k_G12 operators.
Definition: vrr_11_r12kg12_11.h:35
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
Definition: stdarray.h:18
Definition: prefactors.h:159
GenOper is a single operator described by descriptor Descr.
Definition: oper.h:152
Generic integral over a two-body operator with one bfs for each particle in bra and ket...
Definition: integral_11_11.h:32
RRImpl must inherit GenericRecurrenceRelation<RRImpl>
Definition: generic_rr.h:48
QuantumNumbersA<T,N> is a set of N quantum numbers of type T implemented in terms of a C-style array...
Definition: quanta.h:198
std::string to_string(const T &x)
Converts x to its string representation.
Definition: entity.h:71
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:91
Set of basis functions.
Definition: bfset.h:41
bool zero() const
norm() == 0
Definition: bfset.h:153
Helps GenericRecurrenceRelation to work around the compiler problem with make_child.
Definition: generic_rr.h:148