20 #ifndef _libint2_src_bin_libint_itr11twoprep11_h_ 21 #define _libint2_src_bin_libint_itr11twoprep11_h_ 31 #include <twoprep_11_11.h> 34 #include <prefactors.h> 36 #include <default_params.h> 49 template <
template <
typename,
typename,
typename>
class ERI,
class BFSet,
int part, FunctionPosition where>
57 typedef ERI<BFSet,TwoPRep,mType> TargetType;
58 typedef TargetType ChildType;
69 static SafePtr<ThisType> Instance(
const SafePtr<TargetType>&,
unsigned int dir = 0);
82 if (boost::is_same<BasisFunctionType,CGShell>::value)
88 unsigned int num_children()
const {
return children_.size(); };
91 SafePtr<DGVertex>
rr_target()
const {
return static_pointer_cast<
DGVertex,TargetType>(target_); }
93 SafePtr<DGVertex>
rr_child(
unsigned int i)
const {
return static_pointer_cast<
DGVertex,ChildType>(children_.at(i)); }
108 static const unsigned int max_nchildren_ = 4;
110 SafePtr<TargetType> target_;
111 std::vector< SafePtr<ChildType> > children_;
112 const SafePtr<ChildType>& make_child(
const BFSet& A,
const BFSet& B,
const BFSet& C,
const BFSet& D,
unsigned int m) {
113 const SafePtr<ChildType>& i = ChildType::Instance(A,B,C,D,m);
114 children_.push_back(i);
115 return *(children_.end()-1);
118 std::string generate_label()
const 120 typedef typename TargetType::AuxIndexType
mType;
121 static SafePtr<mType> aux0(
new mType(0u));
125 os <<
"ITR Part" << part <<
" " <<
to_string(where) << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
129 #if LIBINT_ENABLE_GENERIC_CODE 130 bool has_generic(
const SafePtr<CompilationParameters>& cparams)
const;
133 std::string generic_header()
const;
135 std::string generic_instance(
const SafePtr<CodeContext>& context,
const SafePtr<CodeSymbols>& args)
const;
139 template <
template <
typename,
typename,
typename>
class ERI,
class F,
int part, FunctionPosition where>
140 SafePtr< ITR_11_TwoPRep_11<ERI,F,part,where> >
144 SafePtr<ThisType> this_ptr(
new ThisType(Tint,dir));
146 if (this_ptr->num_children() != 0) {
147 this_ptr->template register_with_rrstack<ThisType>();
150 return SafePtr<ThisType>();
153 template <
template <
typename,
typename,
typename>
class ERI,
class F,
int part, FunctionPosition where>
156 target_(Tint), dir_(dir)
164 if (a.contracted() ||
177 const bool deriv = dA.
zero() ==
false ||
178 dB.
zero() ==
false ||
179 dC.
zero() ==
false ||
185 children_.reserve(max_nchildren_);
188 const unsigned int m = Tint->aux()->elem(0);
189 const F& _1 = unit<F>(dir);
192 if (part == 0 && where == InBra) {
193 F a(Tint->bra(0,0) - _1);
200 if (not (b.zero() && d.zero()))
return;
202 const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
203 if (is_simple()) { expr_ = Vector(
"TwoPRepITR_pfac0_0_0")[dir] * ABCD_m; nflops_+=1; }
205 const F& cp1 = c + _1;
206 const SafePtr<ChildType>& ABCp1D_m = make_child(a,b,cp1,d,m);
207 if (is_simple()) { expr_ -= Scalar(
"eoz") * ABCp1D_m; nflops_+=2; }
209 const F& am1 = a - _1;
211 const SafePtr<ChildType>& Am1BCD_m = make_child(am1,b,c,d,m);
212 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2z") * Am1BCD_m; nflops_+=3; }
214 const F& cm1 = c - _1;
216 const SafePtr<ChildType>& ABCm1D_m = make_child(a,b,cm1,d,m);
217 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2z") * ABCm1D_m; nflops_+=3; }
222 if (part == 1 && where == InBra) {
225 F c(Tint->bra(1,0) - _1);
230 if (not (b.zero() && d.zero()))
return;
232 const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
233 if (is_simple()) { expr_ = Vector(
"TwoPRepITR_pfac0_1_0")[dir] * ABCD_m; nflops_+=1; }
235 const F& ap1 = a + _1;
236 const SafePtr<ChildType>& Ap1BCD_m = make_child(ap1,b,c,d,m);
237 if (is_simple()) { expr_ -= Scalar(
"zoe") * Ap1BCD_m; nflops_+=2; }
239 const F& cm1 = c - _1;
241 const SafePtr<ChildType>& ABCm1D_m = make_child(a,b,cm1,d,m);
242 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar(
"oo2e") * ABCm1D_m; nflops_+=3; }
244 const F& am1 = a - _1;
246 const SafePtr<ChildType>& Am1BCD_m = make_child(am1,b,c,d,m);
247 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar(
"oo2e") * Am1BCD_m; nflops_+=3; }
252 if (part == 0 && where == InKet) {
254 F b(Tint->ket(0,0) - _1);
260 if (not (a.zero() && c.zero()))
return;
262 const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
263 if (is_simple()) { expr_ = Vector(
"TwoPRepITR_pfac0_0_1")[dir] * ABCD_m; nflops_+=1; }
265 const F& dp1 = d + _1;
266 const SafePtr<ChildType>& ABCDp1_m = make_child(a,b,c,dp1,m);
267 if (is_simple()) { expr_ -= Scalar(
"eoz") * ABCDp1_m; nflops_+=2; }
269 const F& bm1 = b - _1;
271 const SafePtr<ChildType>& ABm1CD_m = make_child(a,bm1,c,d,m);
272 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2z") * ABm1CD_m; nflops_+=3; }
274 const F& dm1 = d - _1;
276 const SafePtr<ChildType>& ABCDm1_m = make_child(a,b,c,dm1,m);
277 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2z") * ABCDm1_m; nflops_+=3; }
282 if (part == 1 && where == InKet) {
286 F d(Tint->ket(1,0) - _1);
290 if (not (a.zero() && c.zero()))
return;
292 const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
293 if (is_simple()) { expr_ = Vector(
"TwoPRepITR_pfac0_1_1")[dir] * ABCD_m; nflops_+=1; }
295 const F& bp1 = b + _1;
296 const SafePtr<ChildType>& ABp1CD_m = make_child(a,bp1,c,d,m);
297 if (is_simple()) { expr_ -= Scalar(
"zoe") * ABp1CD_m; nflops_+=2; }
299 const F& dm1 = d - _1;
301 const SafePtr<ChildType>& ABCDm1_m = make_child(a,b,c,dm1,m);
302 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar(
"oo2e") * ABCDm1_m; nflops_+=3; }
304 const F& bm1 = b - _1;
306 const SafePtr<ChildType>& ABm1CD_m = make_child(a,bm1,c,d,m);
307 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar(
"oo2e") * ABm1CD_m; nflops_+=3; }
314 #if LIBINT_ENABLE_GENERIC_CODE 315 template <
template <
typename,
typename,
typename>
class ERI,
class F,
int part, FunctionPosition where>
322 F sh_a(target_->bra(0,0));
323 F sh_b(target_->ket(0,0));
324 F sh_c(target_->bra(1,0));
325 F sh_d(target_->ket(1,0));
326 const unsigned int max_opt_am = cparams->max_am_opt();
329 if (sh_b.zero() && sh_d.zero() &&
330 (sh_a.norm() > std::max(2*max_opt_am,1u) ||
331 sh_c.norm() > std::max(2*max_opt_am,1u)
333 (sh_a.norm() > 1u && sh_c.norm() > 1u)
335 const bool ITR_xs_xs_Part1_implemented =
false;
336 if (part == 1)
return ITR_xs_xs_Part1_implemented;
339 if (sh_a.zero() && sh_c.zero() &&
340 (sh_b.norm() > std::max(2*max_opt_am,1u) ||
341 sh_d.norm() > std::max(2*max_opt_am,1u)
343 (sh_b.norm() > 1u && sh_d.norm() > 1u)
345 const bool ITR_sx_sx_implemented =
false;
346 return ITR_sx_sx_implemented;
351 template <
template <
typename,
typename,
typename>
class ERI,
class F,
int part, FunctionPosition where>
355 F sh_a(target_->bra(0,0));
356 F sh_b(target_->ket(0,0));
357 F sh_c(target_->bra(1,0));
358 F sh_d(target_->ket(1,0));
359 const bool xsxs = sh_b.zero() && sh_d.zero();
360 const bool sxsx = sh_a.zero() && sh_c.zero();
366 const bool deriv = dA.
zero() ==
false ||
367 dB.
zero() ==
false ||
368 dC.
zero() ==
false ||
370 assert(deriv ==
false);
372 if (deriv ==
false) {
373 return std::string(
"ITR_xs_xs.h");
376 if (xsxs)
return std::string(
"ITR_xs_xs_deriv.h");
377 if (sxsx)
return std::string(
"ITR_sx_sx_deriv.h");
382 template <
template <
typename,
typename,
typename>
class ERI,
class F,
int part, FunctionPosition where>
386 std::ostringstream oss;
387 F sh_a(target_->bra(0,0));
388 F sh_b(target_->ket(0,0));
389 F sh_c(target_->bra(1,0));
390 F sh_d(target_->ket(1,0));
391 const bool xsxs = sh_b.zero() && sh_d.zero();
392 const bool sxsx = sh_a.zero() && sh_c.zero();
398 const bool deriv = dA.
zero() ==
false ||
399 dB.
zero() ==
false ||
400 dC.
zero() ==
false ||
402 assert(deriv ==
false);
404 oss <<
"using namespace libint2;" << endl;
406 if (deriv ==
false) {
408 oss <<
"libint2::ITR_xs_xs<" << part <<
"," << sh_a.norm() <<
"," << sh_c.norm() <<
",true,";
411 oss <<
"libint2::ITR_xs_xs<" << part <<
"," << sh_b.norm() <<
"," << sh_d.norm() <<
",false,";
413 oss << ((context->cparams()->max_vector_length() == 1) ?
"false" :
"true");
414 oss <<
">::compute(inteval";
416 const unsigned int nargs = args->n();
417 for(
unsigned int a=0; a<nargs; a++) {
418 oss <<
"," << args->symbol(a);
425 #endif // #if !LIBINT_ENABLE_GENERIC_CODE TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition: bfset.h:879
int partindex_direction() const
overrides RecurrenceRelation::partindex_direction()
Definition: itr_11_twoprep_11.h:73
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
Definition: stdarray.h:18
Definition: prefactors.h:159
bool is_simple() const
Implementation of RecurrenceRelation::is_simple()
Definition: itr_11_twoprep_11.h:95
This is a vertex of a Directed Graph (DG)
Definition: dgvertex.h:42
SafePtr< DGVertex > rr_child(unsigned int i) const
Implementation of RecurrenceRelation::rr_child()
Definition: itr_11_twoprep_11.h:93
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
AlgebraicOperator is an algebraic operator that acts on objects of type T.
Definition: algebra.h:47
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:91
RecurrenceRelation::ExprType ExprType
The type of expressions in which RecurrenceRelations result.
Definition: itr_11_twoprep_11.h:60
ITR (Interelectron Transfer Relation) for 2-e ERI.
Definition: itr_11_twoprep_11.h:50
RecurrenceRelation describes all recurrence relations.
Definition: rr.h:100
Set of basis functions.
Definition: bfset.h:41
SafePtr< DGVertex > rr_target() const
Implementation of RecurrenceRelation::rr_target()
Definition: itr_11_twoprep_11.h:91
static bool directional()
Default directionality.
Definition: itr_11_twoprep_11.h:81
bool zero() const
norm() == 0
Definition: bfset.h:153