LIBINT  2.6.0
vrr_11_r12kg12_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_vrr11r12kg1211_h_
22 #define _libint2_src_bin_libint_vrr11r12kg1211_h_
23 
24 #include <generic_rr.h>
25 #include <r12kg12_11_11.h>
26 
27 using namespace std;
28 
29 namespace libint2 {
30 
35  template <class BFSet, int part, FunctionPosition where>
36  class VRR_11_R12kG12_11 : public GenericRecurrenceRelation< VRR_11_R12kG12_11<BFSet,part,where>,
37  BFSet,
38  GenIntegralSet_11_11<BFSet,R12kG12,mType> >
39  {
40  public:
42  typedef BFSet BasisFunctionType;
46  static const unsigned int max_nchildren = 8;
47 
48  using ParentType::Instance;
49 
51  static bool directional() { return ParentType::default_directional(); }
52 
53  private:
54  using ParentType::RecurrenceRelation::expr_;
55  using ParentType::RecurrenceRelation::nflops_;
56  using ParentType::target_;
57  using ParentType::is_simple;
58 
60  VRR_11_R12kG12_11(const SafePtr<TargetType>&, unsigned int dir);
61 
62  static std::string descr() { return "VRR"; }
66  std::string generate_label() const
67  {
68  typedef typename TargetType::AuxIndexType mType;
69  static SafePtr<mType> aux0(new mType(0u));
70  ostringstream os;
71  os << descr() << "P" << part << to_string(where)
72  << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
73  return os.str();
74  }
75 
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;
83  #endif
84  };
85 
86  template <class F, int part, FunctionPosition where>
87  VRR_11_R12kG12_11<F,part,where>::VRR_11_R12kG12_11(const SafePtr< TargetType >& Tint,
88  unsigned int dir) :
89  ParentType(Tint,dir)
90  {
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);
98 
99  // does not work for contracted functions
100  {
101  F a(Tint->bra(0,0));
102  F b(Tint->ket(0,0));
103  F c(Tint->bra(1,0));
104  F d(Tint->ket(1,0));
105  if (a.contracted() ||
106  b.contracted() ||
107  c.contracted() ||
108  d.contracted() ||
109  Tint->oper()->descr().contracted())
110  return;
111  }
112 
113  // does not work for derivative integrals (yet or ever)
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 ||
121  dD.zero() == false;
122  assert(deriv == false);
123 
124 
125  typedef TargetType ChildType;
126  ChildFactory<ThisType,ChildType> factory(this);
127 
128  // if K is -1, the recurrence relation looks exactly as it would for ERI
129  // thus generate the same code, and remember to use appropriate prefactors
130  if (K == -1) {
131  // Build on A
132  if (part == 0 && where == InBra) {
133  F a(Tint->bra(0,0) - _1);
134  if (!exists(a)) return;
135  F b(Tint->ket(0,0));
136  F c(Tint->bra(1,0));
137  F d(Tint->ket(1,0));
138 
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; }
142 
143  const F& am1 = a - _1;
144  if (exists(am1)) {
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; }
148  }
149  const F& bm1 = b - _1;
150  if (exists(bm1)) {
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; }
154  }
155  const F& cm1 = c - _1;
156  if (exists(cm1)) {
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; }
159  }
160  const F& dm1 = d - _1;
161  if (exists(dm1)) {
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; }
164  }
165  return;
166  }
167  // Build on A
168  if (part == 0 && where == InKet) {
169  F a(Tint->bra(0,0));
170  F b(Tint->ket(0,0) - _1);
171  if (!exists(b)) return;
172  F c(Tint->bra(1,0));
173  F d(Tint->ket(1,0));
174 
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; }
178 
179  const F& am1 = a - _1;
180  if (exists(am1)) {
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; }
184  }
185  const F& bm1 = b - _1;
186  if (exists(bm1)) {
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; }
190  }
191  const F& cm1 = c - _1;
192  if (exists(cm1)) {
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; }
195  }
196  const F& dm1 = d - _1;
197  if (exists(dm1)) {
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; }
200  }
201  return;
202  }
203  // Build on C
204  if (part == 1 && where == InBra) {
205  F a(Tint->bra(0,0));
206  F b(Tint->ket(0,0));
207  F c(Tint->bra(1,0) - _1);
208  if (!exists(c)) return;
209  F d(Tint->ket(1,0));
210 
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; }
214 
215  const F& cm1 = c - _1;
216  if (exists(cm1)) {
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; }
220  }
221  const F& dm1 = d - _1;
222  if (exists(dm1)) {
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; }
226  }
227  const F& am1 = a - _1;
228  if (exists(am1)) {
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; }
231  }
232  const F& bm1 = b - _1;
233  if (exists(bm1)) {
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; }
236  }
237  return;
238  }
239  // Build on D
240  if (part == 1 && where == InKet) {
241  F a(Tint->bra(0,0));
242  F b(Tint->ket(0,0));
243  F c(Tint->bra(1,0));
244  F d(Tint->ket(1,0) - _1);
245  if (!exists(d)) return;
246 
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; }
250 
251  const F& cm1 = c - _1;
252  if (exists(cm1)) {
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; }
256  }
257  const F& dm1 = d - _1;
258  if (exists(dm1)) {
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; }
262  }
263  const F& am1 = a - _1;
264  if (exists(am1)) {
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; }
267  }
268  const F& bm1 = b - _1;
269  if (exists(bm1)) {
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; }
272  }
273  return;
274  }
275  return;
276  } // K == -1?
277  else {
278  // K != -1, the auxiliary quantum number is not used
279  if (m != 0)
280  throw std::logic_error("VRR_11_R12kG12_11<I,F,K,part,where>::children_and_expr_Kge0() -- nonzero auxiliary quantum detected.");
281 
282  // can build (a0|c0) or (0b|0d)
283  bool xsxs = false;
284  bool sxsx = false;
285  {
286  F b(Tint->ket(0,0) - _1);
287  F d(Tint->ket(1,0) - _1);
288  if (!exists(b) && !exists(d))
289  xsxs = true;
290  }
291  {
292  F a(Tint->bra(0,0) - _1);
293  F c(Tint->bra(1,0) - _1);
294  if (!exists(a) && !exists(c))
295  sxsx = true;
296  }
297  // can't handle the general case
298  if (!xsxs && !sxsx)
299  return;
300  // can't handle (ss|ss) case
301  if (xsxs && sxsx)
302  return;
303 
304  if (xsxs) {
305  // Build on A
306  if (part == 0 && where == InBra) {
307  F a(Tint->bra(0,0) - _1);
308  if (!exists(a)) return;
309  F b(Tint->ket(0,0));
310  F c(Tint->bra(1,0));
311  F d(Tint->ket(1,0));
312 
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;
316  if (exists(am1)) {
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; }
319  }
320  const F& cm1 = c - _1;
321  if (exists(cm1)) {
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; }
324  }
325  if (K != 0) {
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; }
332  }
333  return;
334  }
335  // Build on C
336  if (part == 1 && where == InBra) {
337  F a(Tint->bra(0,0));
338  F b(Tint->ket(0,0));
339  F c(Tint->bra(1,0) - _1);
340  if (!exists(c)) return;
341  F d(Tint->ket(1,0));
342 
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;
346  if (exists(cm1)) {
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; }
349  }
350  const F& am1 = a - _1;
351  if (exists(am1)) {
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; }
354  }
355  if (K != 0) {
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; }
362  }
363  return;
364  }
365  } // end of a0c0 case
366 
367  if (sxsx) {
368  // Build on B
369  if (part == 0 && where == InKet) {
370  F a(Tint->bra(0,0));
371  F b(Tint->ket(0,0) - _1);
372  if (!exists(b)) return;
373  F c(Tint->bra(1,0));
374  F d(Tint->ket(1,0));
375 
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;
379  if (exists(bm1)) {
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; }
382  }
383  const F& dm1 = d - _1;
384  if (exists(dm1)) {
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; }
387  }
388  if (K != 0) {
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; }
395  }
396  return;
397  }
398  // Build on D
399  if (part == 1 && where == InKet) {
400  F a(Tint->bra(0,0));
401  F b(Tint->ket(0,0));
402  F c(Tint->bra(1,0));
403  F d(Tint->ket(1,0) - _1);
404  if (!exists(d)) return;
405 
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;
409  if (exists(dm1)) {
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; }
412  }
413  const F& bm1 = b - _1;
414  if (exists(bm1)) {
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; }
417  }
418  if (K != 0) {
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; }
425  }
426  return;
427  }
428  } // end of 0b0d case
429 
430  return;
431  } // K >= 0
432  }
433 
434  #if LIBINT_ENABLE_GENERIC_CODE
435  template <class F, int part, FunctionPosition where>
436  bool
437  VRR_11_R12kG12_11<F,part,where>::has_generic(const SafePtr<CompilationParameters>& cparams) const
438  {
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();
444  // generic code works for a0c0 and 0a0c classes where am(a) > 1 and am(c) > 1
445  // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
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)
450  ) &&
451  (sh_a.norm() > 1u && sh_c.norm() > 1u)
452  )
453  return true;
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)
458  ) &&
459  (sh_b.norm() > 1u && sh_d.norm() > 1u)
460  )
461  return true;
462  return false;
463  }
464 
465  template <class F, int part, FunctionPosition where>
466  std::string
467  VRR_11_R12kG12_11<F,part,where>::generic_header() const
468  {
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();
476  if (K == -1) {
477  if (xsxs)
478  return std::string("OSVRR_xs_xs.h");
479  if (sxsx)
480  return std::string("OSVRR_sx_sx.h");
481  }
482  else {
483  return std::string("VRR_r12kg12_xs_xs.h");
484  }
485  assert(false); // unreachable
486  }
487 
488  template <class F, int part, FunctionPosition where>
489  std::string
490  VRR_11_R12kG12_11<F,part,where>::generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const
491  {
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();
500 
501  bool ahlrichs_simplification = false;
502  bool unit_s = false;
503  if (xsxs) {
504  unit_s = sh_b.is_unit();
505  }
506  else { // (sxsx)
507  unit_s = sh_a.is_unit();
508  }
509 
510  oss << "using namespace libint2;" << endl;
511  if (K == -1) {
512  if (xsxs) {
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");
516  }
517  if (sxsx) {
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");
521  }
522  }
523  else {
524  if (xsxs) {
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");
527  }
528  if (sxsx) {
529  // NOTE that using same function, the only difference is in the RR prefactors
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");
532  }
533  }
534  oss << ">::compute(inteval";
535 
536  // if K == -1 follow the same logic as for standard VRR
537  if (K == -1) {
538  oss << "," << args->symbol(0); // target
539  if (not ahlrichs_simplification && unit_s) // purely to avoid having a 4-term generic RR, reuse 5-term with dummy argument
540  oss << ",0"; // src0-> 0x0
541  const unsigned int nargs = args->n();
542  for(unsigned int a=1; a<nargs; a++) {
543  oss << "," << args->symbol(a);
544  }
545  }
546  else { // K != -1
547  const unsigned int nargs = args->n();
548  for(unsigned int a=0; a<nargs; a++) {
549  oss << "," << args->symbol(a);
550  }
551  }
552  // if K == 0 add dummy arguments so that the same generic function can be used for all K>=0 cases
553  if (K == 0) {
554  for(unsigned int a=0; a<3; ++a) {
555  oss << ",(LIBINT2_REALTYPE*)0";
556  }
557  }
558 
559  oss << ");";
560 
561  return oss.str();
562  }
563  #endif // #if !LIBINT_ENABLE_GENERIC_CODE
564 
565 };
566 
567 #endif
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