LIBINT  2.1.0-stable
itr_11_twoprep_11.h
1 /*
2  * This file is a part of Libint.
3  * Copyright (C) 2004-2014 Edward F. Valeev
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 2 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see http://www.gnu.org/licenses/.
17  *
18  */
19 
20 #ifndef _libint2_src_bin_libint_itr11twoprep11_h_
21 #define _libint2_src_bin_libint_itr11twoprep11_h_
22 
23 #include <iostream>
24 #include <sstream>
25 #include <string>
26 #include <vector>
27 #include <stdexcept>
28 #include <assert.h>
29 #include <dgvertex.h>
30 #include <rr.h>
31 #include <twoprep_11_11.h>
32 #include <algebra.h>
33 #include <flop.h>
34 #include <prefactors.h>
35 #include <context.h>
36 #include <default_params.h>
37 #include <util.h>
38 
39 using namespace std;
40 
41 
42 namespace libint2 {
43 
49  template <template <typename,typename,typename> class ERI, class BFSet, int part, FunctionPosition where>
51  {
52 
53  public:
55  typedef BFSet BasisFunctionType;
57  typedef ERI<BFSet,TwoPRep,mType> TargetType;
58  typedef TargetType ChildType;
61 
69  static SafePtr<ThisType> Instance(const SafePtr<TargetType>&, unsigned int dir = 0);
70  virtual ~ITR_11_TwoPRep_11() { assert(part == 0 || part == 1); }
71 
73  int partindex_direction() const { return part == 0 ? +1 // transfer from 0 to 1
74  : -1; // transfer from 1 to 0
75  }
76 
78 
81  static bool directional() {
82  if (boost::is_same<BasisFunctionType,CGShell>::value)
83  return false;
84  return true;
85  }
86 
87 #if 1
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)); }
95  bool is_simple() const {
97  }
98 #endif
99 
100  private:
106  ITR_11_TwoPRep_11(const SafePtr<TargetType>&, unsigned int dir);
107 
108  static const unsigned int max_nchildren_ = 4;
109  unsigned int dir_;
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);
116  }
117 
118  std::string generate_label() const
119  {
120  typedef typename TargetType::AuxIndexType mType;
121  static SafePtr<mType> aux0(new mType(0u));
122  ostringstream os;
123  // ITR recurrence relation code is independent of m (it never appears anywhere in equations), hence
124  // to avoid generating identical code make sure that the (unique) label does not contain m
125  os << "ITR Part" << part << " " << to_string(where) << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
126  return os.str();
127  }
128 
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;
136 #endif
137  };
138 
139  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
140  SafePtr< ITR_11_TwoPRep_11<ERI,F,part,where> >
141  ITR_11_TwoPRep_11<ERI,F,part,where>::Instance(const SafePtr<TargetType>& Tint,
142  unsigned int dir)
143  {
144  SafePtr<ThisType> this_ptr(new ThisType(Tint,dir));
145  // Do post-construction duties
146  if (this_ptr->num_children() != 0) {
147  this_ptr->template register_with_rrstack<ThisType>();
148  return this_ptr;
149  }
150  return SafePtr<ThisType>();
151  }
152 
153  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
154  ITR_11_TwoPRep_11<ERI,F,part,where>::ITR_11_TwoPRep_11(const SafePtr<TargetType>& Tint,
155  unsigned int dir) :
156  target_(Tint), dir_(dir)
157  {
158  // only works for primitive integrals
159  {
160  F a(Tint->bra(0,0));
161  F b(Tint->ket(0,0));
162  F c(Tint->bra(1,0));
163  F d(Tint->ket(1,0));
164  if (a.contracted() ||
165  b.contracted() ||
166  c.contracted() ||
167  d.contracted())
168  return;
169  }
170 
171  // not yet implemented for derivative integrals
172  {
173  const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
174  const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
175  const OriginDerivative<3u> dC = Tint->bra(1,0).deriv();
176  const OriginDerivative<3u> dD = Tint->ket(1,0).deriv();
177  const bool deriv = dA.zero() == false ||
178  dB.zero() == false ||
179  dC.zero() == false ||
180  dD.zero() == false;
181  if (deriv)
182  return;
183  }
184 
185  children_.reserve(max_nchildren_);
186  using namespace libint2::algebra;
187  using namespace libint2::prefactor;
188  const unsigned int m = Tint->aux()->elem(0);
189  const F& _1 = unit<F>(dir);
190 
191  // Build on A
192  if (part == 0 && where == InBra) {
193  F a(Tint->bra(0,0) - _1);
194  if (!exists(a)) return;
195  F b(Tint->ket(0,0));
196  F c(Tint->bra(1,0));
197  F d(Tint->ket(1,0));
198 
199  // must be (A0|C0)
200  if (not (b.zero() && d.zero())) return;
201 
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; }
204 
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; }
208 
209  const F& am1 = a - _1;
210  if (exists(am1)) {
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; }
213  }
214  const F& cm1 = c - _1;
215  if (exists(cm1)) {
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; }
218  }
219  return;
220  }
221  // Build on C
222  if (part == 1 && where == InBra) {
223  F a(Tint->bra(0,0));
224  F b(Tint->ket(0,0));
225  F c(Tint->bra(1,0) - _1);
226  if (!exists(c)) return;
227  F d(Tint->ket(1,0));
228 
229  // must be (A0|C0)
230  if (not (b.zero() && d.zero())) return;
231 
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; }
234 
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; }
238 
239  const F& cm1 = c - _1;
240  if (exists(cm1)) {
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; }
243  }
244  const F& am1 = a - _1;
245  if (exists(am1)) {
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; }
248  }
249  return;
250  }
251  // Build on B
252  if (part == 0 && where == InKet) {
253  F a(Tint->bra(0,0));
254  F b(Tint->ket(0,0) - _1);
255  if (!exists(b)) return;
256  F c(Tint->bra(1,0));
257  F d(Tint->ket(1,0));
258 
259  // must be (0B|0D)
260  if (not (a.zero() && c.zero())) return;
261 
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; }
264 
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; }
268 
269  const F& bm1 = b - _1;
270  if (exists(bm1)) {
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; }
273  }
274  const F& dm1 = d - _1;
275  if (exists(dm1)) {
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; }
278  }
279  return;
280  }
281  // Build on D
282  if (part == 1 && where == InKet) {
283  F a(Tint->bra(0,0));
284  F b(Tint->ket(0,0));
285  F c(Tint->bra(1,0));
286  F d(Tint->ket(1,0) - _1);
287  if (!exists(d)) return;
288 
289  // must be (0B|0D)
290  if (not (a.zero() && c.zero())) return;
291 
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; }
294 
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; }
298 
299  const F& dm1 = d - _1;
300  if (exists(dm1)) {
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; }
303  }
304  const F& bm1 = b - _1;
305  if (exists(bm1)) {
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; }
308  }
309  return;
310  }
311 
312  }
313 
314 #if LIBINT_ENABLE_GENERIC_CODE
315  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
316  bool
317  ITR_11_TwoPRep_11<ERI,F,part,where>::has_generic(const SafePtr<CompilationParameters>& cparams) const
318  {
320  return false;
321 
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();
327  // generic code works for a0c0 of 0a0c classes where am(a) > 1 and am(c) > 1
328  // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
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)
332  ) &&
333  (sh_a.norm() > 1u && sh_c.norm() > 1u)
334  ) {
335  const bool ITR_xs_xs_Part1_implemented = false; // only Part0 is implemented
336  if (part == 1) return ITR_xs_xs_Part1_implemented;
337  else return true;
338  }
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)
342  ) &&
343  (sh_b.norm() > 1u && sh_d.norm() > 1u)
344  ) {
345  const bool ITR_sx_sx_implemented = false; // only ITR_xs_xs is implemented
346  return ITR_sx_sx_implemented;
347  }
348  return false;
349  }
350 
351  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
352  std::string
354  {
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();
361 
362  const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
363  const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
364  const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
365  const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
366  const bool deriv = dA.zero() == false ||
367  dB.zero() == false ||
368  dC.zero() == false ||
369  dD.zero() == false;
370  assert(deriv == false);
371 
372  if (deriv == false) {
373  return std::string("ITR_xs_xs.h");
374  }
375  else {
376  if (xsxs) return std::string("ITR_xs_xs_deriv.h");
377  if (sxsx) return std::string("ITR_sx_sx_deriv.h");
378  }
379  abort(); // unreachable
380  }
381 
382  template <template <typename,typename,typename> class ERI, class F, int part, FunctionPosition where>
383  std::string
384  ITR_11_TwoPRep_11<ERI,F,part,where>::generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const
385  {
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();
393 
394  const OriginDerivative<3u> dA = target_->bra(0,0).deriv();
395  const OriginDerivative<3u> dB = target_->ket(0,0).deriv();
396  const OriginDerivative<3u> dC = target_->bra(1,0).deriv();
397  const OriginDerivative<3u> dD = target_->ket(1,0).deriv();
398  const bool deriv = dA.zero() == false ||
399  dB.zero() == false ||
400  dC.zero() == false ||
401  dD.zero() == false;
402  assert(deriv == false);
403 
404  oss << "using namespace libint2;" << endl;
405 
406  if (deriv == false) { // for regular integrals I know exactly how many prerequisites I need
407  if(xsxs) {
408  oss << "libint2::ITR_xs_xs<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",true,";
409  }
410  if (sxsx) {
411  oss << "libint2::ITR_xs_xs<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",false,";
412  }
413  oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
414  oss << ">::compute(inteval";
415 
416  const unsigned int nargs = args->n();
417  for(unsigned int a=0; a<nargs; a++) {
418  oss << "," << args->symbol(a);
419  }
420  oss << ");";
421  }
422 
423  return oss.str();
424  }
425 #endif // #if !LIBINT_ENABLE_GENERIC_CODE
426 
427 };
428 
429 #endif
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
Definition: algebra.h:32
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