LIBINT  2.6.0
itr_11_twoprep_11.h
1 /*
2  * Copyright (C) 2004-2019 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 
40 using namespace std;
41 
42 
43 namespace libint2 {
44 
50  template <template <typename,typename,typename> class ERI, class BFSet, int part, FunctionPosition where>
52  {
53 
54  public:
56  typedef BFSet BasisFunctionType;
58  typedef ERI<BFSet,TwoPRep,mType> TargetType;
59  typedef TargetType ChildType;
62 
70  static SafePtr<ThisType> Instance(const SafePtr<TargetType>&, unsigned int dir = 0);
71  virtual ~ITR_11_TwoPRep_11() { assert(part == 0 || part == 1); }
72 
74  int partindex_direction() const { return part == 0 ? +1 // transfer from 0 to 1
75  : -1; // transfer from 1 to 0
76  }
77 
79 
82  static bool directional() {
83  if (boost::is_same<BasisFunctionType,CGShell>::value)
84  return false;
85  return true;
86  }
87 
88 #if 1
89  unsigned int num_children() const { return children_.size(); };
92  SafePtr<DGVertex> rr_target() const { return static_pointer_cast<DGVertex,TargetType>(target_); }
94  SafePtr<DGVertex> rr_child(unsigned int i) const { return static_pointer_cast<DGVertex,ChildType>(children_.at(i)); }
96  bool is_simple() const {
98  }
99 #endif
100 
101  private:
107  ITR_11_TwoPRep_11(const SafePtr<TargetType>&, unsigned int dir);
108 
109  static const unsigned int max_nchildren_ = 4;
110  unsigned int dir_;
111  SafePtr<TargetType> target_;
112  std::vector< SafePtr<ChildType> > children_;
113  const SafePtr<ChildType>& make_child(const BFSet& A, const BFSet& B, const BFSet& C, const BFSet& D, unsigned int m) {
114  const SafePtr<ChildType>& i = ChildType::Instance(A,B,C,D,m);
115  children_.push_back(i);
116  return *(children_.end()-1);
117  }
118 
119  std::string generate_label() const
120  {
121  typedef typename TargetType::AuxIndexType mType;
122  static SafePtr<mType> aux0(new mType(0u));
123  ostringstream os;
124  // ITR recurrence relation code is independent of m (it never appears anywhere in equations), hence
125  // to avoid generating identical code make sure that the (unique) label does not contain m
126  os << "ITR Part" << part << " " << to_string(where) << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
127  return os.str();
128  }
129 
130 #if LIBINT_ENABLE_GENERIC_CODE
131  bool has_generic(const SafePtr<CompilationParameters>& cparams) const;
134  std::string generic_header() const;
136  std::string generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const;
137 #endif
138  };
139 
140  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
141  SafePtr< ITR_11_TwoPRep_11<ERI,F,part,where> >
142  ITR_11_TwoPRep_11<ERI,F,part,where>::Instance(const SafePtr<TargetType>& Tint,
143  unsigned int dir)
144  {
145  SafePtr<ThisType> this_ptr(new ThisType(Tint,dir));
146  // Do post-construction duties
147  if (this_ptr->num_children() != 0) {
148  this_ptr->template register_with_rrstack<ThisType>();
149  return this_ptr;
150  }
151  return SafePtr<ThisType>();
152  }
153 
154  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
155  ITR_11_TwoPRep_11<ERI,F,part,where>::ITR_11_TwoPRep_11(const SafePtr<TargetType>& Tint,
156  unsigned int dir) :
157  target_(Tint), dir_(dir)
158  {
159  // only works for primitive integrals
160  {
161  F a(Tint->bra(0,0));
162  F b(Tint->ket(0,0));
163  F c(Tint->bra(1,0));
164  F d(Tint->ket(1,0));
165  if (a.contracted() ||
166  b.contracted() ||
167  c.contracted() ||
168  d.contracted())
169  return;
170  }
171 
172  // not yet implemented for derivative integrals
173  {
174  const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
175  const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
176  const OriginDerivative<3u> dC = Tint->bra(1,0).deriv();
177  const OriginDerivative<3u> dD = Tint->ket(1,0).deriv();
178  const bool deriv = dA.zero() == false ||
179  dB.zero() == false ||
180  dC.zero() == false ||
181  dD.zero() == false;
182  if (deriv)
183  return;
184  }
185 
186  children_.reserve(max_nchildren_);
187  using namespace libint2::algebra;
188  using namespace libint2::prefactor;
189  const unsigned int m = Tint->aux()->elem(0);
190  const F& _1 = unit<F>(dir);
191 
192  // Build on A
193  if (part == 0 && where == InBra) {
194  F a(Tint->bra(0,0) - _1);
195  if (!exists(a)) return;
196  F b(Tint->ket(0,0));
197  F c(Tint->bra(1,0));
198  F d(Tint->ket(1,0));
199 
200  // must be (A0|C0)
201  if (not (b.zero() && d.zero())) return;
202 
203  const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
204  if (is_simple()) { expr_ = Vector("TwoPRepITR_pfac0_0_0")[dir] * ABCD_m; nflops_+=1; }
205 
206  const F& cp1 = c + _1;
207  const SafePtr<ChildType>& ABCp1D_m = make_child(a,b,cp1,d,m);
208  if (is_simple()) { expr_ -= Scalar("eoz") * ABCp1D_m; nflops_+=2; }
209 
210  const F& am1 = a - _1;
211  if (exists(am1)) {
212  const SafePtr<ChildType>& Am1BCD_m = make_child(am1,b,c,d,m);
213  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * Am1BCD_m; nflops_+=3; }
214  }
215  const F& cm1 = c - _1;
216  if (exists(cm1)) {
217  const SafePtr<ChildType>& ABCm1D_m = make_child(a,b,cm1,d,m);
218  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2z") * ABCm1D_m; nflops_+=3; }
219  }
220  return;
221  }
222  // Build on C
223  if (part == 1 && where == InBra) {
224  F a(Tint->bra(0,0));
225  F b(Tint->ket(0,0));
226  F c(Tint->bra(1,0) - _1);
227  if (!exists(c)) return;
228  F d(Tint->ket(1,0));
229 
230  // must be (A0|C0)
231  if (not (b.zero() && d.zero())) return;
232 
233  const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
234  if (is_simple()) { expr_ = Vector("TwoPRepITR_pfac0_1_0")[dir] * ABCD_m; nflops_+=1; }
235 
236  const F& ap1 = a + _1;
237  const SafePtr<ChildType>& Ap1BCD_m = make_child(ap1,b,c,d,m);
238  if (is_simple()) { expr_ -= Scalar("zoe") * Ap1BCD_m; nflops_+=2; }
239 
240  const F& cm1 = c - _1;
241  if (exists(cm1)) {
242  const SafePtr<ChildType>& ABCm1D_m = make_child(a,b,cm1,d,m);
243  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2e") * ABCm1D_m; nflops_+=3; }
244  }
245  const F& am1 = a - _1;
246  if (exists(am1)) {
247  const SafePtr<ChildType>& Am1BCD_m = make_child(am1,b,c,d,m);
248  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2e") * Am1BCD_m; nflops_+=3; }
249  }
250  return;
251  }
252  // Build on B
253  if (part == 0 && where == InKet) {
254  F a(Tint->bra(0,0));
255  F b(Tint->ket(0,0) - _1);
256  if (!exists(b)) return;
257  F c(Tint->bra(1,0));
258  F d(Tint->ket(1,0));
259 
260  // must be (0B|0D)
261  if (not (a.zero() && c.zero())) return;
262 
263  const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
264  if (is_simple()) { expr_ = Vector("TwoPRepITR_pfac0_0_1")[dir] * ABCD_m; nflops_+=1; }
265 
266  const F& dp1 = d + _1;
267  const SafePtr<ChildType>& ABCDp1_m = make_child(a,b,c,dp1,m);
268  if (is_simple()) { expr_ -= Scalar("eoz") * ABCDp1_m; nflops_+=2; }
269 
270  const F& bm1 = b - _1;
271  if (exists(bm1)) {
272  const SafePtr<ChildType>& ABm1CD_m = make_child(a,bm1,c,d,m);
273  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * ABm1CD_m; nflops_+=3; }
274  }
275  const F& dm1 = d - _1;
276  if (exists(dm1)) {
277  const SafePtr<ChildType>& ABCDm1_m = make_child(a,b,c,dm1,m);
278  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2z") * ABCDm1_m; nflops_+=3; }
279  }
280  return;
281  }
282  // Build on D
283  if (part == 1 && where == InKet) {
284  F a(Tint->bra(0,0));
285  F b(Tint->ket(0,0));
286  F c(Tint->bra(1,0));
287  F d(Tint->ket(1,0) - _1);
288  if (!exists(d)) return;
289 
290  // must be (0B|0D)
291  if (not (a.zero() && c.zero())) return;
292 
293  const SafePtr<ChildType>& ABCD_m = make_child(a,b,c,d,m);
294  if (is_simple()) { expr_ = Vector("TwoPRepITR_pfac0_1_1")[dir] * ABCD_m; nflops_+=1; }
295 
296  const F& bp1 = b + _1;
297  const SafePtr<ChildType>& ABp1CD_m = make_child(a,bp1,c,d,m);
298  if (is_simple()) { expr_ -= Scalar("zoe") * ABp1CD_m; nflops_+=2; }
299 
300  const F& dm1 = d - _1;
301  if (exists(dm1)) {
302  const SafePtr<ChildType>& ABCDm1_m = make_child(a,b,c,dm1,m);
303  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2e") * ABCDm1_m; nflops_+=3; }
304  }
305  const F& bm1 = b - _1;
306  if (exists(bm1)) {
307  const SafePtr<ChildType>& ABm1CD_m = make_child(a,bm1,c,d,m);
308  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2e") * ABm1CD_m; nflops_+=3; }
309  }
310  return;
311  }
312 
313  }
314 
315 #if LIBINT_ENABLE_GENERIC_CODE
316  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
317  bool
318  ITR_11_TwoPRep_11<ERI,F,part,where>::has_generic(const SafePtr<CompilationParameters>& cparams) const
319  {
320  if (TrivialBFSet<F>::result)
321  return false;
322 
323  F sh_a(target_->bra(0,0));
324  F sh_b(target_->ket(0,0));
325  F sh_c(target_->bra(1,0));
326  F sh_d(target_->ket(1,0));
327  const unsigned int max_opt_am = cparams->max_am_opt();
328  // generic code works for a0c0 of 0a0c classes where am(a) > 1 and am(c) > 1
329  // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
330  if (sh_b.zero() && sh_d.zero() &&
331  (sh_a.norm() > std::max(2*max_opt_am,1u) ||
332  sh_c.norm() > std::max(2*max_opt_am,1u)
333  ) &&
334  (sh_a.norm() > 1u && sh_c.norm() > 1u)
335  ) {
336  const bool ITR_xs_xs_Part1_implemented = false; // only Part0 is implemented
337  if (part == 1) return ITR_xs_xs_Part1_implemented;
338  else return true;
339  }
340  if (sh_a.zero() && sh_c.zero() &&
341  (sh_b.norm() > std::max(2*max_opt_am,1u) ||
342  sh_d.norm() > std::max(2*max_opt_am,1u)
343  ) &&
344  (sh_b.norm() > 1u && sh_d.norm() > 1u)
345  ) {
346  const bool ITR_sx_sx_implemented = false; // only ITR_xs_xs is implemented
347  return ITR_sx_sx_implemented;
348  }
349  return false;
350  }
351 
352  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
353  std::string
354  ITR_11_TwoPRep_11<ERI,F,part,where>::generic_header() const
355  {
356  F sh_a(target_->bra(0,0));
357  F sh_b(target_->ket(0,0));
358  F sh_c(target_->bra(1,0));
359  F sh_d(target_->ket(1,0));
360  const bool xsxs = sh_b.zero() && sh_d.zero();
361  const bool sxsx = sh_a.zero() && sh_c.zero();
362 
363  const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
364  const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
365  const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
366  const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
367  const bool deriv = dA.zero() == false ||
368  dB.zero() == false ||
369  dC.zero() == false ||
370  dD.zero() == false;
371  assert(deriv == false);
372 
373  if (deriv == false) {
374  return std::string("ITR_xs_xs.h");
375  }
376  else {
377  if (xsxs) return std::string("ITR_xs_xs_deriv.h");
378  if (sxsx) return std::string("ITR_sx_sx_deriv.h");
379  }
380  abort(); // unreachable
381  }
382 
383  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
384  std::string
385  ITR_11_TwoPRep_11<ERI,F,part,where>::generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const
386  {
387  std::ostringstream oss;
388  F sh_a(target_->bra(0,0));
389  F sh_b(target_->ket(0,0));
390  F sh_c(target_->bra(1,0));
391  F sh_d(target_->ket(1,0));
392  const bool xsxs = sh_b.zero() && sh_d.zero();
393  const bool sxsx = sh_a.zero() && sh_c.zero();
394 
395  const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
396  const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
397  const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
398  const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
399  const bool deriv = dA.zero() == false ||
400  dB.zero() == false ||
401  dC.zero() == false ||
402  dD.zero() == false;
403  assert(deriv == false);
404 
405  oss << "using namespace libint2;" << endl;
406 
407  if (deriv == false) { // for regular integrals I know exactly how many prerequisites I need
408  if(xsxs) {
409  oss << "libint2::ITR_xs_xs<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",true,";
410  }
411  if (sxsx) {
412  oss << "libint2::ITR_xs_xs<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",false,";
413  }
414  oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
415  oss << ">::compute(inteval";
416 
417  const unsigned int nargs = args->n();
418  for(unsigned int a=0; a<nargs; a++) {
419  oss << "," << args->symbol(a);
420  }
421  oss << ");";
422  }
423 
424  return oss.str();
425  }
426 #endif // #if !LIBINT_ENABLE_GENERIC_CODE
427 
428 };
429 
430 #endif
SafePtr< DGVertex > rr_child(unsigned int i) const
Implementation of RecurrenceRelation::rr_child()
Definition: itr_11_twoprep_11.h:94
TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition: bfset.h:892
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
bool is_simple() const
Implementation of RecurrenceRelation::is_simple()
Definition: itr_11_twoprep_11.h:96
AlgebraicOperator is an algebraic operator that acts on objects of type T.
Definition: algebra.h:48
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
RecurrenceRelation::ExprType ExprType
The type of expressions in which RecurrenceRelations result.
Definition: itr_11_twoprep_11.h:61
ITR (Interelectron Transfer Relation) for 2-e ERI.
Definition: itr_11_twoprep_11.h:51
RecurrenceRelation describes all recurrence relations.
Definition: rr.h:101
Set of basis functions.
Definition: bfset.h:42
static bool directional()
Default directionality.
Definition: itr_11_twoprep_11.h:82
int partindex_direction() const
overrides RecurrenceRelation::partindex_direction()
Definition: itr_11_twoprep_11.h:74
SafePtr< DGVertex > rr_target() const
Implementation of RecurrenceRelation::rr_target()
Definition: itr_11_twoprep_11.h:92