21 #ifndef _libint2_src_bin_libint_vrr11r12kg1211_h_ 22 #define _libint2_src_bin_libint_vrr11r12kg1211_h_ 24 #include <generic_rr.h> 25 #include <r12kg12_11_11.h> 35 template <
class BFSet,
int part, FunctionPosition where>
38 GenIntegralSet_11_11<BFSet,R12kG12,mType> >
46 static const unsigned int max_nchildren = 8;
48 using ParentType::Instance;
51 static bool directional() {
return ParentType::default_directional(); }
54 using ParentType::RecurrenceRelation::expr_;
55 using ParentType::RecurrenceRelation::nflops_;
56 using ParentType::target_;
57 using ParentType::is_simple;
62 static std::string descr() {
return "VRR"; }
66 std::string generate_label()
const 68 typedef typename TargetType::AuxIndexType mType;
69 static SafePtr<mType> aux0(
new mType(0u));
71 os << descr() <<
"P" << part << to_string(where)
72 << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
76 #if LIBINT_ENABLE_GENERIC_CODE 77 bool has_generic(
const SafePtr<CompilationParameters>& cparams)
const;
80 std::string generic_header()
const;
82 std::string generic_instance(
const SafePtr<CodeContext>& context,
const SafePtr<CodeSymbols>& args)
const;
86 template <
class F,
int part, FunctionPosition where>
87 VRR_11_R12kG12_11<F,part,where>::VRR_11_R12kG12_11(
const SafePtr< TargetType >& Tint,
91 using namespace libint2::algebra;
92 using namespace libint2::prefactor;
93 const int K = Tint->oper()->descr().K();
94 typedef R12_k_G12_Descr R12kG12Descr;
95 const R12kG12 oK((R12kG12Descr(K)));
96 const unsigned int m = Tint->aux()->elem(0);
97 const F _1 = unit<F>(dir);
105 if (a.contracted() ||
109 Tint->oper()->descr().contracted())
114 const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
115 const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
116 const OriginDerivative<3u> dC = Tint->bra(1,0).deriv();
117 const OriginDerivative<3u> dD = Tint->ket(1,0).deriv();
118 const bool deriv = dA.zero() ==
false ||
119 dB.zero() ==
false ||
120 dC.zero() ==
false ||
122 assert(deriv ==
false);
125 typedef TargetType ChildType;
126 ChildFactory<ThisType,ChildType> factory(
this);
132 if (part == 0 && where == InBra) {
133 F a(Tint->bra(0,0) - _1);
139 auto ABCD_m = factory.make_child(a,b,c,d,m,oK);
140 auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1,oK);
141 if (
is_simple()) { expr_ = Vector(
"PA")[dir] * ABCD_m + Vector(
"WP")[dir] * ABCD_mp1; nflops_+=3; }
143 const F& am1 = a - _1;
145 auto Am1BCD_m = factory.make_child(am1,b,c,d,m,oK);
146 auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1,oK);
147 if (
is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2z") * (Am1BCD_m - Scalar(
"roz") * Am1BCD_mp1); nflops_+=5; }
149 const F& bm1 = b - _1;
151 auto ABm1CD_m = factory.make_child(a,bm1,c,d,m,oK);
152 auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1,oK);
153 if (
is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2z") * (ABm1CD_m - Scalar(
"roz") * ABm1CD_mp1); nflops_+=5; }
155 const F& cm1 = c - _1;
157 auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1,oK);
158 if (
is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2ze") * ABCm1D_mp1; nflops_+=3; }
160 const F& dm1 = d - _1;
162 auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1,oK);
163 if (
is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2ze") * ABCDm1_mp1; nflops_+=3; }
168 if (part == 0 && where == InKet) {
170 F b(Tint->ket(0,0) - _1);
175 auto ABCD_m = factory.make_child(a,b,c,d,m,oK);
176 auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1,oK);
177 if (
is_simple()) { expr_ = Vector(
"PB")[dir] * ABCD_m + Vector(
"WP")[dir] * ABCD_mp1; nflops_+=3; }
179 const F& am1 = a - _1;
181 auto Am1BCD_m = factory.make_child(am1,b,c,d,m,oK);
182 auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1,oK);
183 if (
is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2z") * (Am1BCD_m - Scalar(
"roz") * Am1BCD_mp1); nflops_+=5; }
185 const F& bm1 = b - _1;
187 auto ABm1CD_m = factory.make_child(a,bm1,c,d,m,oK);
188 auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1,oK);
189 if (
is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2z") * (ABm1CD_m - Scalar(
"roz") * ABm1CD_mp1); nflops_+=5; }
191 const F& cm1 = c - _1;
193 auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1,oK);
194 if (
is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2ze") * ABCm1D_mp1; nflops_+=3; }
196 const F& dm1 = d - _1;
198 auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1,oK);
199 if (
is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2ze") * ABCDm1_mp1; nflops_+=3; }
204 if (part == 1 && where == InBra) {
207 F c(Tint->bra(1,0) - _1);
211 auto ABCD_m = factory.make_child(a,b,c,d,m,oK);
212 auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1,oK);
213 if (
is_simple()) { expr_ = Vector(
"QC")[dir] * ABCD_m + Vector(
"WQ")[dir] * ABCD_mp1; nflops_+=3; }
215 const F& cm1 = c - _1;
217 auto ABCm1D_m = factory.make_child(a,b,cm1,d,m,oK);
218 auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1,oK);
219 if (
is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2e") * (ABCm1D_m - Scalar(
"roe") * ABCm1D_mp1); nflops_+=5; }
221 const F& dm1 = d - _1;
223 auto ABCDm1_m = factory.make_child(a,b,c,dm1,m,oK);
224 auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1,oK);
225 if (
is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2e") * (ABCDm1_m - Scalar(
"roe") * ABCDm1_mp1); nflops_+=5; }
227 const F& am1 = a - _1;
229 auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1,oK);
230 if (
is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2ze") * Am1BCD_mp1; nflops_+=3; }
232 const F& bm1 = b - _1;
234 auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1,oK);
235 if (
is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2ze") * ABm1CD_mp1; nflops_+=3; }
240 if (part == 1 && where == InKet) {
244 F d(Tint->ket(1,0) - _1);
247 auto ABCD_m = factory.make_child(a,b,c,d,m,oK);
248 auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1,oK);
249 if (
is_simple()) { expr_ = Vector(
"QD")[dir] * ABCD_m + Vector(
"WQ")[dir] * ABCD_mp1; nflops_+=3; }
251 const F& cm1 = c - _1;
253 auto ABCm1D_m = factory.make_child(a,b,cm1,d,m,oK);
254 auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1,oK);
255 if (
is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2e") * (ABCm1D_m - Scalar(
"roe") * ABCm1D_mp1); nflops_+=5; }
257 const F& dm1 = d - _1;
259 auto ABCDm1_m = factory.make_child(a,b,c,dm1,m,oK);
260 auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1,oK);
261 if (
is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2e") * (ABCDm1_m - Scalar(
"roe") * ABCDm1_mp1); nflops_+=5; }
263 const F& am1 = a - _1;
265 auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1,oK);
266 if (
is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2ze") * Am1BCD_mp1; nflops_+=3; }
268 const F& bm1 = b - _1;
270 auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1,oK);
271 if (
is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2ze") * ABm1CD_mp1; nflops_+=3; }
280 throw std::logic_error(
"VRR_11_R12kG12_11<I,F,K,part,where>::children_and_expr_Kge0() -- nonzero auxiliary quantum detected.");
286 F b(Tint->ket(0,0) - _1);
287 F d(Tint->ket(1,0) - _1);
292 F a(Tint->bra(0,0) - _1);
293 F c(Tint->bra(1,0) - _1);
306 if (part == 0 && where == InBra) {
307 F a(Tint->bra(0,0) - _1);
313 auto ABCD_K = factory.make_child(a,b,c,d,0u,oK);
314 if (
is_simple()) { expr_ = Vector(
"R12kG12_pfac0_0")[dir] * ABCD_K; nflops_+=1; }
315 const F& am1 = a - _1;
317 auto Am1BCD_K = factory.make_child(am1,b,c,d,0u,oK);
318 if (
is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"R12kG12_pfac1_0") * Am1BCD_K; nflops_+=3; }
320 const F& cm1 = c - _1;
322 auto ABCm1D_K = factory.make_child(a,b,cm1,d,0u,oK);
323 if (
is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"R12kG12_pfac2") * ABCm1D_K; nflops_+=3; }
326 const R12kG12 oKm2((R12kG12Descr(K-2)));
327 auto Ap1BCD_Km2 = factory.make_child(a+_1,b,c,d,0u,oKm2);
328 auto ABCp1D_Km2 = factory.make_child(a,b,c+_1,d,0u,oKm2);
329 auto ABCD_Km2 = factory.make_child(a,b,c,d,0u,oKm2);
330 if (
is_simple()) { expr_ += Scalar((
double)K) * Scalar(
"R12kG12_pfac3_0")
331 * (Ap1BCD_Km2 - ABCp1D_Km2 + Vector(
"R12kG12_pfac4_0")[dir] * ABCD_Km2); nflops_+=6; }
336 if (part == 1 && where == InBra) {
339 F c(Tint->bra(1,0) - _1);
343 auto ABCD_K = factory.make_child(a,b,c,d,0u,oK);
344 if (
is_simple()) { expr_ = Vector(
"R12kG12_pfac0_1")[dir] * ABCD_K; nflops_+=1; }
345 const F& cm1 = c - _1;
347 auto ABCm1D_K = factory.make_child(a,b,cm1,d,0u,oK);
348 if (
is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"R12kG12_pfac1_1") * ABCm1D_K; nflops_+=3; }
350 const F& am1 = a - _1;
352 auto Am1BCD_K = factory.make_child(am1,b,c,d,0u,oK);
353 if (
is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"R12kG12_pfac2") * Am1BCD_K; nflops_+=3; }
356 const R12kG12 oKm2((R12kG12Descr(K-2)));
357 auto ABCp1D_Km2 = factory.make_child(a,b,c+_1,d,0u,oKm2);
358 auto Ap1BCD_Km2 = factory.make_child(a+_1,b,c,d,0u,oKm2);
359 auto ABCD_Km2 = factory.make_child(a,b,c,d,0u,oKm2);
360 if (
is_simple()) { expr_ += Scalar((
double)K) * Scalar(
"R12kG12_pfac3_1")
361 * (ABCp1D_Km2 - Ap1BCD_Km2 + Vector(
"R12kG12_pfac4_1")[dir] * ABCD_Km2); nflops_+=6; }
369 if (part == 0 && where == InKet) {
371 F b(Tint->ket(0,0) - _1);
376 auto ABCD_K = factory.make_child(a,b,c,d,0u,oK);
377 if (
is_simple()) { expr_ = Vector(
"R12kG12_pfac0_0")[dir] * ABCD_K; nflops_+=1; }
378 const F& bm1 = b - _1;
380 auto ABm1CD_K = factory.make_child(a,bm1,c,d,0u,oK);
381 if (
is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"R12kG12_pfac1_0") * ABm1CD_K; nflops_+=3; }
383 const F& dm1 = d - _1;
385 auto ABCDm1_K = factory.make_child(a,b,c,dm1,0u,oK);
386 if (
is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"R12kG12_pfac2") * ABCDm1_K; nflops_+=3; }
389 const R12kG12 oKm2((R12kG12Descr(K-2)));
390 auto ABp1CD_Km2 = factory.make_child(a,b+_1,c,d,0u,oKm2);
391 auto ABCDp1_Km2 = factory.make_child(a,b,c,d+_1,0u,oKm2);
392 auto ABCD_Km2 = factory.make_child(a,b,c,d,0u,oKm2);
393 if (
is_simple()) { expr_ += Scalar((
double)K) * Scalar(
"R12kG12_pfac3_0")
394 * (ABp1CD_Km2 - ABCDp1_Km2 + Vector(
"R12kG12_pfac4_0")[dir] * ABCD_Km2); nflops_+=6; }
399 if (part == 1 && where == InKet) {
403 F d(Tint->ket(1,0) - _1);
406 auto ABCD_K = factory.make_child(a,b,c,d,0u,oK);
407 if (
is_simple()) { expr_ = Vector(
"R12kG12_pfac0_1")[dir] * ABCD_K; nflops_+=1; }
408 const F& dm1 = d - _1;
410 auto ABCDm1_K = factory.make_child(a,b,c,dm1,0u,oK);
411 if (
is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"R12kG12_pfac1_1") * ABCDm1_K; nflops_+=3; }
413 const F& bm1 = b - _1;
415 auto ABm1CD_K = factory.make_child(a,bm1,c,d,0u,oK);
416 if (
is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"R12kG12_pfac2") * ABm1CD_K; nflops_+=3; }
419 const R12kG12 oKm2((R12kG12Descr(K-2)));
420 auto ABCDp1_Km2 = factory.make_child(a,b,c,d+_1,0u,oKm2);
421 auto ABp1CD_Km2 = factory.make_child(a,b+_1,c,d,0u,oKm2);
422 auto ABCD_Km2 = factory.make_child(a,b,c,d,0u,oKm2);
423 if (
is_simple()) { expr_ += Scalar((
double)K) * Scalar(
"R12kG12_pfac3_1")
424 * (ABCDp1_Km2 - ABp1CD_Km2 + Vector(
"R12kG12_pfac4_1")[dir] * ABCD_Km2); nflops_+=6; }
434 #if LIBINT_ENABLE_GENERIC_CODE 435 template <
class F,
int part, FunctionPosition where>
437 VRR_11_R12kG12_11<F,part,where>::has_generic(
const SafePtr<CompilationParameters>& cparams)
const 439 F sh_a(target_->bra(0,0));
440 F sh_b(target_->ket(0,0));
441 F sh_c(target_->bra(1,0));
442 F sh_d(target_->ket(1,0));
443 const unsigned int max_opt_am = cparams->max_am_opt();
446 if (!TrivialBFSet<F>::result &&
447 sh_b.zero() && sh_d.zero() &&
448 (sh_a.norm() > std::max(2*max_opt_am,1u) ||
449 sh_c.norm() > std::max(2*max_opt_am,1u)
451 (sh_a.norm() > 1u && sh_c.norm() > 1u)
454 if (!TrivialBFSet<F>::result &&
455 sh_a.zero() && sh_c.zero() &&
456 (sh_b.norm() > std::max(2*max_opt_am,1u) ||
457 sh_d.norm() > std::max(2*max_opt_am,1u)
459 (sh_b.norm() > 1u && sh_d.norm() > 1u)
465 template <
class F,
int part, FunctionPosition where>
467 VRR_11_R12kG12_11<F,part,where>::generic_header()
const 469 F sh_a(target_->bra(0,0));
470 F sh_b(target_->ket(0,0));
471 F sh_c(target_->bra(1,0));
472 F sh_d(target_->ket(1,0));
473 const bool xsxs = sh_b.zero() && sh_d.zero();
474 const bool sxsx = sh_a.zero() && sh_c.zero();
475 const int K = target_->oper()->descr().K();
478 return std::string(
"OSVRR_xs_xs.h");
480 return std::string(
"OSVRR_sx_sx.h");
483 return std::string(
"VRR_r12kg12_xs_xs.h");
488 template <
class F,
int part, FunctionPosition where>
490 VRR_11_R12kG12_11<F,part,where>::generic_instance(
const SafePtr<CodeContext>& context,
const SafePtr<CodeSymbols>& args)
const 492 const int K = target_->oper()->descr().K();
493 std::ostringstream oss;
494 F sh_a(target_->bra(0,0));
495 F sh_b(target_->ket(0,0));
496 F sh_c(target_->bra(1,0));
497 F sh_d(target_->ket(1,0));
498 const bool xsxs = sh_b.zero() && sh_d.zero();
499 const bool sxsx = sh_a.zero() && sh_c.zero();
501 bool ahlrichs_simplification =
false;
504 unit_s = sh_b.is_unit();
507 unit_s = sh_a.is_unit();
510 oss <<
"using namespace libint2;" << endl;
513 oss <<
"libint2::OSVRR_xs_xs<" << part <<
"," << sh_a.norm() <<
"," << sh_c.norm() <<
",";
514 if (not ahlrichs_simplification) oss << (unit_s ?
"true," :
"false,");
515 oss << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
518 oss <<
"libint2::OSVRR_sx_sx<" << part <<
"," << sh_b.norm() <<
"," << sh_d.norm() <<
",";
519 if (not ahlrichs_simplification) oss << (unit_s ?
"true," :
"false,");
520 oss << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
525 oss <<
"libint2::VRR_r12kg12_xs_xs<" << part <<
"," << sh_a.norm() <<
"," << sh_c.norm() <<
",";
526 oss << K <<
"," << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
530 oss <<
"libint2::VRR_r12kg12_xs_xs<" << part <<
"," << sh_b.norm() <<
"," << sh_d.norm() <<
",";
531 oss << K <<
"," << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
534 oss <<
">::compute(inteval";
538 oss <<
"," << args->symbol(0);
539 if (not ahlrichs_simplification && unit_s)
541 const unsigned int nargs = args->n();
542 for(
unsigned int a=1; a<nargs; a++) {
543 oss <<
"," << args->symbol(a);
547 const unsigned int nargs = args->n();
548 for(
unsigned int a=0; a<nargs; a++) {
549 oss <<
"," << args->symbol(a);
554 for(
unsigned int a=0; a<3; ++a) {
555 oss <<
",(LIBINT2_REALTYPE*)0";
563 #endif // #if !LIBINT_ENABLE_GENERIC_CODE static bool directional()
Default directionality.
Definition: vrr_11_r12kg12_11.h:51
VRR Recurrence Relation for 2-e integrals of the R12_k_G12 operators.
Definition: vrr_11_r12kg12_11.h:36
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
bool is_simple() const
Implementation of RecurrenceRelation::is_simple()
Definition: generic_rr.h:81
Generic integral over a two-body operator with one bfs for each particle in bra and ket.
Definition: integral_11_11.h:33
RRImpl must inherit GenericRecurrenceRelation<RRImpl>
Definition: generic_rr.h:49
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
Set of basis functions.
Definition: bfset.h:42