LIBINT 2.7.2
itr_11_twoprep_11.h
1/*
2 * Copyright (C) 2004-2021 Edward F. Valeev
3 *
4 * This file is part of Libint.
5 *
6 * Libint is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * Libint is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with Libint. If not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21#ifndef _libint2_src_bin_libint_itr11twoprep11_h_
22#define _libint2_src_bin_libint_itr11twoprep11_h_
23
24#include <iostream>
25#include <sstream>
26#include <string>
27#include <vector>
28#include <stdexcept>
29#include <cassert>
30#include <dgvertex.h>
31#include <rr.h>
32#include <twoprep_11_11.h>
33#include <algebra.h>
34#include <flop.h>
35#include <prefactors.h>
36#include <context.h>
37#include <default_params.h>
38#include <util.h>
39
40namespace libint2 {
41
47 template <template <typename,typename,typename> class ERI, class BFSet, int part, FunctionPosition where>
48 class ITR_11_TwoPRep_11 : public RecurrenceRelation
49 {
50
51 public:
52 typedef RecurrenceRelation ParentType;
53 typedef BFSet BasisFunctionType;
54 typedef ITR_11_TwoPRep_11 ThisType;
55 typedef ERI<BFSet,TwoPRep,mType> TargetType;
56 typedef TargetType ChildType;
58 typedef RecurrenceRelation::ExprType ExprType;
59
67 static SafePtr<ThisType> Instance(const SafePtr<TargetType>&, unsigned int dir = 0);
68 virtual ~ITR_11_TwoPRep_11() { assert(part == 0 || part == 1); }
69
71 int partindex_direction() const { return part == 0 ? +1 // transfer from 0 to 1
72 : -1; // transfer from 1 to 0
73 }
74
76
79 static bool directional() {
80 if (boost::is_same<BasisFunctionType,CGShell>::value)
81 return false;
82 return true;
83 }
84
85#if 1
87 unsigned int num_children() const override { return children_.size(); }
89 SafePtr<DGVertex> rr_target() const override { return static_pointer_cast<DGVertex,TargetType>(target_); }
91 SafePtr<DGVertex> rr_child(unsigned int i) const override { return static_pointer_cast<DGVertex,ChildType>(children_.at(i)); }
93 bool is_simple() const override {
94 return TrivialBFSet<BFSet>::result;
95 }
96#endif
97
98 private:
104 ITR_11_TwoPRep_11(const SafePtr<TargetType>&, unsigned int dir);
105
106 static const unsigned int max_nchildren_ = 4;
107 unsigned int dir_;
108 SafePtr<TargetType> target_;
109 std::vector< SafePtr<ChildType> > children_;
110 const SafePtr<ChildType>& make_child(const BFSet& A, const BFSet& B, const BFSet& C, const BFSet& D, unsigned int m) {
111 const SafePtr<ChildType>& i = ChildType::Instance(A,B,C,D,m);
112 children_.push_back(i);
113 return *(children_.end()-1);
114 }
115
116 std::string generate_label() const override
117 {
118 typedef typename TargetType::AuxIndexType mType;
119 static SafePtr<mType> aux0(new mType(0u));
120 std::ostringstream os;
121 // ITR recurrence relation code is independent of m (it never appears anywhere in equations), hence
122 // to avoid generating identical code make sure that the (unique) label does not contain m
123 os << "ITR Part" << part << " " << to_string(where) << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
124 return os.str();
125 }
126
127#if LIBINT_ENABLE_GENERIC_CODE
129 bool has_generic(const SafePtr<CompilationParameters>& cparams) const override;
131 std::string generic_header() const override;
133 std::string generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const override;
134#endif
135 };
136
137 template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
138 SafePtr< ITR_11_TwoPRep_11<ERI,F,part,where> >
139 ITR_11_TwoPRep_11<ERI,F,part,where>::Instance(const SafePtr<TargetType>& Tint,
140 unsigned int dir)
141 {
142 SafePtr<ThisType> this_ptr(new ThisType(Tint,dir));
143 // Do post-construction duties
144 if (this_ptr->num_children() != 0) {
145 this_ptr->template register_with_rrstack<ThisType>();
146 return this_ptr;
147 }
148 return SafePtr<ThisType>();
149 }
150
151 template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
152 ITR_11_TwoPRep_11<ERI,F,part,where>::ITR_11_TwoPRep_11(const SafePtr<TargetType>& Tint,
153 unsigned int dir) :
154 target_(Tint), dir_(dir)
155 {
156 // only works for primitive integrals
157 {
158 F a(Tint->bra(0,0));
159 F b(Tint->ket(0,0));
160 F c(Tint->bra(1,0));
161 F d(Tint->ket(1,0));
162 if (a.contracted() ||
163 b.contracted() ||
164 c.contracted() ||
165 d.contracted())
166 return;
167 }
168
169 // not yet implemented for derivative integrals
170 {
171 const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
172 const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
173 const OriginDerivative<3u> dC = Tint->bra(1,0).deriv();
174 const OriginDerivative<3u> dD = Tint->ket(1,0).deriv();
175 const bool deriv = !dA.zero() ||
176 !dB.zero() ||
177 !dC.zero() ||
178 !dD.zero();
179 if (deriv)
180 return;
181 }
182
183 children_.reserve(max_nchildren_);
184 using namespace libint2::algebra;
185 using namespace libint2::prefactor;
186 const unsigned int m = Tint->aux()->elem(0);
187 const F& _1 = unit<F>(dir);
188
189 // Build on A
190 if (part == 0 && where == InBra) {
191 F a(Tint->bra(0,0) - _1);
192 if (!exists(a)) return;
193 F b(Tint->ket(0,0));
194 F c(Tint->bra(1,0));
195 F d(Tint->ket(1,0));
196
197 // must be (A0|C0)
198 if (not (b.zero() && d.zero())) return;
199
200 const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
201 if (is_simple()) { expr_ = Vector("TwoPRepITR_pfac0_0_0")[dir] * ABCD_m; nflops_+=1; }
202
203 const F& cp1 = c + _1;
204 const SafePtr<ChildType>& ABCp1D_m = make_child(a,b,cp1,d,m);
205 if (is_simple()) { expr_ -= Scalar("eoz") * ABCp1D_m; nflops_+=2; }
206
207 const F& am1 = a - _1;
208 if (exists(am1)) {
209 const SafePtr<ChildType>& Am1BCD_m = make_child(am1,b,c,d,m);
210 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * Am1BCD_m; nflops_+=3; }
211 }
212 const F& cm1 = c - _1;
213 if (exists(cm1)) {
214 const SafePtr<ChildType>& ABCm1D_m = make_child(a,b,cm1,d,m);
215 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2z") * ABCm1D_m; nflops_+=3; }
216 }
217 return;
218 }
219 // Build on C
220 if (part == 1 && where == InBra) {
221 F a(Tint->bra(0,0));
222 F b(Tint->ket(0,0));
223 F c(Tint->bra(1,0) - _1);
224 if (!exists(c)) return;
225 F d(Tint->ket(1,0));
226
227 // must be (A0|C0)
228 if (not (b.zero() && d.zero())) return;
229
230 const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
231 if (is_simple()) { expr_ = Vector("TwoPRepITR_pfac0_1_0")[dir] * ABCD_m; nflops_+=1; }
232
233 const F& ap1 = a + _1;
234 const SafePtr<ChildType>& Ap1BCD_m = make_child(ap1,b,c,d,m);
235 if (is_simple()) { expr_ -= Scalar("zoe") * Ap1BCD_m; nflops_+=2; }
236
237 const F& cm1 = c - _1;
238 if (exists(cm1)) {
239 const SafePtr<ChildType>& ABCm1D_m = make_child(a,b,cm1,d,m);
240 if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2e") * ABCm1D_m; nflops_+=3; }
241 }
242 const F& am1 = a - _1;
243 if (exists(am1)) {
244 const SafePtr<ChildType>& Am1BCD_m = make_child(am1,b,c,d,m);
245 if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2e") * Am1BCD_m; nflops_+=3; }
246 }
247 return;
248 }
249 // Build on B
250 if (part == 0 && where == InKet) {
251 F a(Tint->bra(0,0));
252 F b(Tint->ket(0,0) - _1);
253 if (!exists(b)) return;
254 F c(Tint->bra(1,0));
255 F d(Tint->ket(1,0));
256
257 // must be (0B|0D)
258 if (not (a.zero() && c.zero())) return;
259
260 const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
261 if (is_simple()) { expr_ = Vector("TwoPRepITR_pfac0_0_1")[dir] * ABCD_m; nflops_+=1; }
262
263 const F& dp1 = d + _1;
264 const SafePtr<ChildType>& ABCDp1_m = make_child(a,b,c,dp1,m);
265 if (is_simple()) { expr_ -= Scalar("eoz") * ABCDp1_m; nflops_+=2; }
266
267 const F& bm1 = b - _1;
268 if (exists(bm1)) {
269 const SafePtr<ChildType>& ABm1CD_m = make_child(a,bm1,c,d,m);
270 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * ABm1CD_m; nflops_+=3; }
271 }
272 const F& dm1 = d - _1;
273 if (exists(dm1)) {
274 const SafePtr<ChildType>& ABCDm1_m = make_child(a,b,c,dm1,m);
275 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2z") * ABCDm1_m; nflops_+=3; }
276 }
277 return;
278 }
279 // Build on D
280 if (part == 1 && where == InKet) {
281 F a(Tint->bra(0,0));
282 F b(Tint->ket(0,0));
283 F c(Tint->bra(1,0));
284 F d(Tint->ket(1,0) - _1);
285 if (!exists(d)) return;
286
287 // must be (0B|0D)
288 if (not (a.zero() && c.zero())) return;
289
290 const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
291 if (is_simple()) { expr_ = Vector("TwoPRepITR_pfac0_1_1")[dir] * ABCD_m; nflops_+=1; }
292
293 const F& bp1 = b + _1;
294 const SafePtr<ChildType>& ABp1CD_m = make_child(a,bp1,c,d,m);
295 if (is_simple()) { expr_ -= Scalar("zoe") * ABp1CD_m; nflops_+=2; }
296
297 const F& dm1 = d - _1;
298 if (exists(dm1)) {
299 const SafePtr<ChildType>& ABCDm1_m = make_child(a,b,c,dm1,m);
300 if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2e") * ABCDm1_m; nflops_+=3; }
301 }
302 const F& bm1 = b - _1;
303 if (exists(bm1)) {
304 const SafePtr<ChildType>& ABm1CD_m = make_child(a,bm1,c,d,m);
305 if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2e") * ABm1CD_m; nflops_+=3; }
306 }
307 return;
308 }
309
310 }
311
312#if LIBINT_ENABLE_GENERIC_CODE
313 template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
314 bool
315 ITR_11_TwoPRep_11<ERI,F,part,where>::has_generic(const SafePtr<CompilationParameters>& cparams) const
316 {
317 if (TrivialBFSet<F>::result)
318 return false;
319
320 F sh_a(target_->bra(0,0));
321 F sh_b(target_->ket(0,0));
322 F sh_c(target_->bra(1,0));
323 F sh_d(target_->ket(1,0));
324 const unsigned int max_opt_am = cparams->max_am_opt();
325 // generic code works for a0c0 of 0a0c classes where am(a) > 1 and am(c) > 1
326 // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
327 if (sh_b.zero() && sh_d.zero() &&
328 (sh_a.norm() > std::max(2*max_opt_am,1u) ||
329 sh_c.norm() > std::max(2*max_opt_am,1u)
330 ) &&
331 (sh_a.norm() > 1u && sh_c.norm() > 1u)
332 ) {
333 const bool ITR_xs_xs_Part1_implemented = false; // only Part0 is implemented
334 if (part == 1) return ITR_xs_xs_Part1_implemented;
335 else return true;
336 }
337 if (sh_a.zero() && sh_c.zero() &&
338 (sh_b.norm() > std::max(2*max_opt_am,1u) ||
339 sh_d.norm() > std::max(2*max_opt_am,1u)
340 ) &&
341 (sh_b.norm() > 1u && sh_d.norm() > 1u)
342 ) {
343 const bool ITR_sx_sx_implemented = false; // only ITR_xs_xs is implemented
344 return ITR_sx_sx_implemented;
345 }
346 return false;
347 }
348
349 template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
350 std::string
351 ITR_11_TwoPRep_11<ERI,F,part,where>::generic_header() const
352 {
353 F sh_a(target_->bra(0,0));
354 F sh_b(target_->ket(0,0));
355 F sh_c(target_->bra(1,0));
356 F sh_d(target_->ket(1,0));
357 const bool xsxs = sh_b.zero() && sh_d.zero();
358 const bool sxsx = sh_a.zero() && sh_c.zero();
359
360 const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
361 const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
362 const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
363 const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
364 const bool deriv = !dA.zero() ||
365 !dB.zero() ||
366 !dC.zero() ||
367 !dD.zero();
368 assert(!deriv);
369
370 if (!deriv) {
371 return std::string("ITR_xs_xs.h");
372 }
373 else {
374 if (xsxs) return std::string("ITR_xs_xs_deriv.h");
375 if (sxsx) return std::string("ITR_sx_sx_deriv.h");
376 }
377 abort(); // unreachable
378 }
379
380 template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
381 std::string
382 ITR_11_TwoPRep_11<ERI,F,part,where>::generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const
383 {
384 std::ostringstream oss;
385 F sh_a(target_->bra(0,0));
386 F sh_b(target_->ket(0,0));
387 F sh_c(target_->bra(1,0));
388 F sh_d(target_->ket(1,0));
389 const bool xsxs = sh_b.zero() && sh_d.zero();
390 const bool sxsx = sh_a.zero() && sh_c.zero();
391
392 const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
393 const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
394 const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
395 const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
396 const bool deriv = !dA.zero() ||
397 !dB.zero() ||
398 !dC.zero() ||
399 !dD.zero();
400 assert(!deriv);
401
402 oss << "using namespace libint2;" << std::endl;
403
404 if (!deriv) { // for regular integrals I know exactly how many prerequisites I need
405 if(xsxs) {
406 oss << "libint2::ITR_xs_xs<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",true,";
407 }
408 if (sxsx) {
409 oss << "libint2::ITR_xs_xs<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",false,";
410 }
411 oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
412 oss << ">::compute(inteval";
413
414 const unsigned int nargs = args->n();
415 for(unsigned int a=0; a<nargs; a++) {
416 oss << "," << args->symbol(a);
417 }
418 oss << ");";
419 }
420
421 return oss.str();
422 }
423#endif // #if !LIBINT_ENABLE_GENERIC_CODE
424
425};
426
427#endif
libint2::AlgebraicOperator
AlgebraicOperator is an algebraic operator that acts on objects of type T.
Definition: algebra.h:50
libint2::BFSet
Set of basis functions.
Definition: bfset.h:42
libint2::ITR_11_TwoPRep_11
ITR (Interelectron Transfer Relation) for 2-e ERI.
Definition: itr_11_twoprep_11.h:49
libint2::ITR_11_TwoPRep_11::is_simple
bool is_simple() const override
Implementation of RecurrenceRelation::is_simple()
Definition: itr_11_twoprep_11.h:93
libint2::ITR_11_TwoPRep_11::num_children
unsigned int num_children() const override
Implementation of RecurrenceRelation::num_children()
Definition: itr_11_twoprep_11.h:87
libint2::ITR_11_TwoPRep_11::partindex_direction
int partindex_direction() const
overrides RecurrenceRelation::partindex_direction()
Definition: itr_11_twoprep_11.h:71
libint2::ITR_11_TwoPRep_11::directional
static bool directional()
Default directionality.
Definition: itr_11_twoprep_11.h:79
libint2::ITR_11_TwoPRep_11::rr_child
SafePtr< DGVertex > rr_child(unsigned int i) const override
Implementation of RecurrenceRelation::rr_child()
Definition: itr_11_twoprep_11.h:91
libint2::ITR_11_TwoPRep_11::ExprType
RecurrenceRelation::ExprType ExprType
The type of expressions in which RecurrenceRelations result.
Definition: itr_11_twoprep_11.h:58
libint2::ITR_11_TwoPRep_11::Instance
static SafePtr< ThisType > Instance(const SafePtr< TargetType > &, unsigned int dir=0)
Use Instance() to obtain an instance of RR.
Definition: itr_11_twoprep_11.h:139
libint2::ITR_11_TwoPRep_11::rr_target
SafePtr< DGVertex > rr_target() const override
Implementation of RecurrenceRelation::rr_target()
Definition: itr_11_twoprep_11.h:89
libint2::RecurrenceRelation
RecurrenceRelation describes all recurrence relations.
Definition: rr.h:99
libint2
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
libint2::exists
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
libint2::mType
DefaultQuantumNumbers< unsignedint, 1 >::Result mType
mType is the type that describes the auxiliary index of standard 2-body repulsion integrals
Definition: quanta.h:408
libint2::to_string
std::string to_string(const T &x)
Converts x to its string representation.
Definition: entity.h:74
libint2::TrivialBFSet
TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition: bfset.h:892
Generated by doxygen 1.9.3