LIBINT  2.1.0-stable
vrr_11_r12kg12_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_vrr11r12kg1211_h_
21 #define _libint2_src_bin_libint_vrr11r12kg1211_h_
22 
23 #include <generic_rr.h>
24 #include <r12kg12_11_11.h>
25 
26 using namespace std;
27 
28 namespace libint2 {
29 
34  template <class BFSet, int part, FunctionPosition where>
35  class VRR_11_R12kG12_11 : public GenericRecurrenceRelation< VRR_11_R12kG12_11<BFSet,part,where>,
36  BFSet,
37  GenIntegralSet_11_11<BFSet,R12kG12,mType> >
38  {
39  public:
41  typedef BFSet BasisFunctionType;
44  friend class GenericRecurrenceRelation<ThisType,BFSet,TargetType>;
45  static const unsigned int max_nchildren = 8;
46 
47  using ParentType::Instance;
48 
50  static bool directional() { return ParentType::default_directional(); }
51 
52  private:
53  using ParentType::RecurrenceRelation::expr_;
54  using ParentType::RecurrenceRelation::nflops_;
55  using ParentType::target_;
56  using ParentType::is_simple;
57 
59  VRR_11_R12kG12_11(const SafePtr<TargetType>&, unsigned int dir);
60 
61  static std::string descr() { return "VRR"; }
65  std::string generate_label() const
66  {
67  typedef typename TargetType::AuxIndexType mType;
68  static SafePtr<mType> aux0(new mType(0u));
69  ostringstream os;
70  os << descr() << "P" << part << to_string(where)
71  << genintegralset_label(target_->bra(),target_->ket(),aux0,target_->oper());
72  return os.str();
73  }
74 
75  #if LIBINT_ENABLE_GENERIC_CODE
76  bool has_generic(const SafePtr<CompilationParameters>& cparams) const;
79  std::string generic_header() const;
81  std::string generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const;
82  #endif
83  };
84 
85  template <class F, int part, FunctionPosition where>
86  VRR_11_R12kG12_11<F,part,where>::VRR_11_R12kG12_11(const SafePtr< TargetType >& Tint,
87  unsigned int dir) :
88  ParentType(Tint,dir)
89  {
90  using namespace libint2::algebra;
91  using namespace libint2::prefactor;
92  const int K = Tint->oper()->descr().K();
93  typedef R12_k_G12_Descr R12kG12Descr;
94  const R12kG12 oK((R12kG12Descr(K)));
95  const unsigned int m = Tint->aux()->elem(0);
96  const F _1 = unit<F>(dir);
97 
98  // does not work for contracted functions
99  {
100  F a(Tint->bra(0,0));
101  F b(Tint->ket(0,0));
102  F c(Tint->bra(1,0));
103  F d(Tint->ket(1,0));
104  if (a.contracted() ||
105  b.contracted() ||
106  c.contracted() ||
107  d.contracted() ||
108  Tint->oper()->descr().contracted())
109  return;
110  }
111 
112  // does not work for derivative integrals (yet or ever)
113  const OriginDerivative<3u> dA = Tint->bra(0,0).deriv();
114  const OriginDerivative<3u> dB = Tint->ket(0,0).deriv();
115  const OriginDerivative<3u> dC = Tint->bra(1,0).deriv();
116  const OriginDerivative<3u> dD = Tint->ket(1,0).deriv();
117  const bool deriv = dA.zero() == false ||
118  dB.zero() == false ||
119  dC.zero() == false ||
120  dD.zero() == false;
121  assert(deriv == false);
122 
123 
124  typedef TargetType ChildType;
125  ChildFactory<ThisType,ChildType> factory(this);
126 
127  // if K is -1, the recurrence relation looks exactly as it would for ERI
128  // thus generate the same code, and remember to use appropriate prefactors
129  if (K == -1) {
130  // Build on A
131  if (part == 0 && where == InBra) {
132  F a(Tint->bra(0,0) - _1);
133  if (!exists(a)) return;
134  F b(Tint->ket(0,0));
135  F c(Tint->bra(1,0));
136  F d(Tint->ket(1,0));
137 
138  auto ABCD_m = factory.make_child(a,b,c,d,m,oK);
139  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1,oK);
140  if (is_simple()) { expr_ = Vector("PA")[dir] * ABCD_m + Vector("WP")[dir] * ABCD_mp1; nflops_+=3; }
141 
142  const F& am1 = a - _1;
143  if (exists(am1)) {
144  auto Am1BCD_m = factory.make_child(am1,b,c,d,m,oK);
145  auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1,oK);
146  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * (Am1BCD_m - Scalar("roz") * Am1BCD_mp1); nflops_+=5; }
147  }
148  const F& bm1 = b - _1;
149  if (exists(bm1)) {
150  auto ABm1CD_m = factory.make_child(a,bm1,c,d,m,oK);
151  auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1,oK);
152  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * (ABm1CD_m - Scalar("roz") * ABm1CD_mp1); nflops_+=5; }
153  }
154  const F& cm1 = c - _1;
155  if (exists(cm1)) {
156  auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1,oK);
157  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2ze") * ABCm1D_mp1; nflops_+=3; }
158  }
159  const F& dm1 = d - _1;
160  if (exists(dm1)) {
161  auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1,oK);
162  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2ze") * ABCDm1_mp1; nflops_+=3; }
163  }
164  return;
165  }
166  // Build on A
167  if (part == 0 && where == InKet) {
168  F a(Tint->bra(0,0));
169  F b(Tint->ket(0,0) - _1);
170  if (!exists(b)) return;
171  F c(Tint->bra(1,0));
172  F d(Tint->ket(1,0));
173 
174  auto ABCD_m = factory.make_child(a,b,c,d,m,oK);
175  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1,oK);
176  if (is_simple()) { expr_ = Vector("PB")[dir] * ABCD_m + Vector("WP")[dir] * ABCD_mp1; nflops_+=3; }
177 
178  const F& am1 = a - _1;
179  if (exists(am1)) {
180  auto Am1BCD_m = factory.make_child(am1,b,c,d,m,oK);
181  auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1,oK);
182  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2z") * (Am1BCD_m - Scalar("roz") * Am1BCD_mp1); nflops_+=5; }
183  }
184  const F& bm1 = b - _1;
185  if (exists(bm1)) {
186  auto ABm1CD_m = factory.make_child(a,bm1,c,d,m,oK);
187  auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1,oK);
188  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2z") * (ABm1CD_m - Scalar("roz") * ABm1CD_mp1); nflops_+=5; }
189  }
190  const F& cm1 = c - _1;
191  if (exists(cm1)) {
192  auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1,oK);
193  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2ze") * ABCm1D_mp1; nflops_+=3; }
194  }
195  const F& dm1 = d - _1;
196  if (exists(dm1)) {
197  auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1,oK);
198  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2ze") * ABCDm1_mp1; nflops_+=3; }
199  }
200  return;
201  }
202  // Build on C
203  if (part == 1 && where == InBra) {
204  F a(Tint->bra(0,0));
205  F b(Tint->ket(0,0));
206  F c(Tint->bra(1,0) - _1);
207  if (!exists(c)) return;
208  F d(Tint->ket(1,0));
209 
210  auto ABCD_m = factory.make_child(a,b,c,d,m,oK);
211  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1,oK);
212  if (is_simple()) { expr_ = Vector("QC")[dir] * ABCD_m + Vector("WQ")[dir] * ABCD_mp1; nflops_+=3; }
213 
214  const F& cm1 = c - _1;
215  if (exists(cm1)) {
216  auto ABCm1D_m = factory.make_child(a,b,cm1,d,m,oK);
217  auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1,oK);
218  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2e") * (ABCm1D_m - Scalar("roe") * ABCm1D_mp1); nflops_+=5; }
219  }
220  const F& dm1 = d - _1;
221  if (exists(dm1)) {
222  auto ABCDm1_m = factory.make_child(a,b,c,dm1,m,oK);
223  auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1,oK);
224  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2e") * (ABCDm1_m - Scalar("roe") * ABCDm1_mp1); nflops_+=5; }
225  }
226  const F& am1 = a - _1;
227  if (exists(am1)) {
228  auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1,oK);
229  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2ze") * Am1BCD_mp1; nflops_+=3; }
230  }
231  const F& bm1 = b - _1;
232  if (exists(bm1)) {
233  auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1,oK);
234  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2ze") * ABm1CD_mp1; nflops_+=3; }
235  }
236  return;
237  }
238  // Build on D
239  if (part == 1 && where == InKet) {
240  F a(Tint->bra(0,0));
241  F b(Tint->ket(0,0));
242  F c(Tint->bra(1,0));
243  F d(Tint->ket(1,0) - _1);
244  if (!exists(d)) return;
245 
246  auto ABCD_m = factory.make_child(a,b,c,d,m,oK);
247  auto ABCD_mp1 = factory.make_child(a,b,c,d,m+1,oK);
248  if (is_simple()) { expr_ = Vector("QD")[dir] * ABCD_m + Vector("WQ")[dir] * ABCD_mp1; nflops_+=3; }
249 
250  const F& cm1 = c - _1;
251  if (exists(cm1)) {
252  auto ABCm1D_m = factory.make_child(a,b,cm1,d,m,oK);
253  auto ABCm1D_mp1 = factory.make_child(a,b,cm1,d,m+1,oK);
254  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("oo2e") * (ABCm1D_m - Scalar("roe") * ABCm1D_mp1); nflops_+=5; }
255  }
256  const F& dm1 = d - _1;
257  if (exists(dm1)) {
258  auto ABCDm1_m = factory.make_child(a,b,c,dm1,m,oK);
259  auto ABCDm1_mp1 = factory.make_child(a,b,c,dm1,m+1,oK);
260  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("oo2e") * (ABCDm1_m - Scalar("roe") * ABCDm1_mp1); nflops_+=5; }
261  }
262  const F& am1 = a - _1;
263  if (exists(am1)) {
264  auto Am1BCD_mp1 = factory.make_child(am1,b,c,d,m+1,oK);
265  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("oo2ze") * Am1BCD_mp1; nflops_+=3; }
266  }
267  const F& bm1 = b - _1;
268  if (exists(bm1)) {
269  auto ABm1CD_mp1 = factory.make_child(a,bm1,c,d,m+1,oK);
270  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("oo2ze") * ABm1CD_mp1; nflops_+=3; }
271  }
272  return;
273  }
274  return;
275  } // K == -1?
276  else {
277  // K != -1, the auxiliary quantum number is not used
278  if (m != 0)
279  throw std::logic_error("VRR_11_R12kG12_11<I,F,K,part,where>::children_and_expr_Kge0() -- nonzero auxiliary quantum detected.");
280 
281  // can build (a0|c0) or (0b|0d)
282  bool xsxs = false;
283  bool sxsx = false;
284  {
285  F b(Tint->ket(0,0) - _1);
286  F d(Tint->ket(1,0) - _1);
287  if (!exists(b) && !exists(d))
288  xsxs = true;
289  }
290  {
291  F a(Tint->bra(0,0) - _1);
292  F c(Tint->bra(1,0) - _1);
293  if (!exists(a) && !exists(c))
294  sxsx = true;
295  }
296  // can't handle the general case
297  if (!xsxs && !sxsx)
298  return;
299  // can't handle (ss|ss) case
300  if (xsxs && sxsx)
301  return;
302 
303  if (xsxs) {
304  // Build on A
305  if (part == 0 && where == InBra) {
306  F a(Tint->bra(0,0) - _1);
307  if (!exists(a)) return;
308  F b(Tint->ket(0,0));
309  F c(Tint->bra(1,0));
310  F d(Tint->ket(1,0));
311 
312  auto ABCD_K = factory.make_child(a,b,c,d,0u,oK);
313  if (is_simple()) { expr_ = Vector("R12kG12_pfac0_0")[dir] * ABCD_K; nflops_+=1; }
314  const F& am1 = a - _1;
315  if (exists(am1)) {
316  auto Am1BCD_K = factory.make_child(am1,b,c,d,0u,oK);
317  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("R12kG12_pfac1_0") * Am1BCD_K; nflops_+=3; }
318  }
319  const F& cm1 = c - _1;
320  if (exists(cm1)) {
321  auto ABCm1D_K = factory.make_child(a,b,cm1,d,0u,oK);
322  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("R12kG12_pfac2") * ABCm1D_K; nflops_+=3; }
323  }
324  if (K != 0) {
325  const R12kG12 oKm2((R12kG12Descr(K-2)));
326  auto Ap1BCD_Km2 = factory.make_child(a+_1,b,c,d,0u,oKm2);
327  auto ABCp1D_Km2 = factory.make_child(a,b,c+_1,d,0u,oKm2);
328  auto ABCD_Km2 = factory.make_child(a,b,c,d,0u,oKm2);
329  if (is_simple()) { expr_ += Scalar((double)K) * Scalar("R12kG12_pfac3_0")
330  * (Ap1BCD_Km2 - ABCp1D_Km2 + Vector("R12kG12_pfac4_0")[dir] * ABCD_Km2); nflops_+=6; }
331  }
332  return;
333  }
334  // Build on C
335  if (part == 1 && where == InBra) {
336  F a(Tint->bra(0,0));
337  F b(Tint->ket(0,0));
338  F c(Tint->bra(1,0) - _1);
339  if (!exists(c)) return;
340  F d(Tint->ket(1,0));
341 
342  auto ABCD_K = factory.make_child(a,b,c,d,0u,oK);
343  if (is_simple()) { expr_ = Vector("R12kG12_pfac0_1")[dir] * ABCD_K; nflops_+=1; }
344  const F& cm1 = c - _1;
345  if (exists(cm1)) {
346  auto ABCm1D_K = factory.make_child(a,b,cm1,d,0u,oK);
347  if (is_simple()) { expr_ += Scalar(c[dir]) * Scalar("R12kG12_pfac1_1") * ABCm1D_K; nflops_+=3; }
348  }
349  const F& am1 = a - _1;
350  if (exists(am1)) {
351  auto Am1BCD_K = factory.make_child(am1,b,c,d,0u,oK);
352  if (is_simple()) { expr_ += Scalar(a[dir]) * Scalar("R12kG12_pfac2") * Am1BCD_K; nflops_+=3; }
353  }
354  if (K != 0) {
355  const R12kG12 oKm2((R12kG12Descr(K-2)));
356  auto ABCp1D_Km2 = factory.make_child(a,b,c+_1,d,0u,oKm2);
357  auto Ap1BCD_Km2 = factory.make_child(a+_1,b,c,d,0u,oKm2);
358  auto ABCD_Km2 = factory.make_child(a,b,c,d,0u,oKm2);
359  if (is_simple()) { expr_ += Scalar((double)K) * Scalar("R12kG12_pfac3_1")
360  * (ABCp1D_Km2 - Ap1BCD_Km2 + Vector("R12kG12_pfac4_1")[dir] * ABCD_Km2); nflops_+=6; }
361  }
362  return;
363  }
364  } // end of a0c0 case
365 
366  if (sxsx) {
367  // Build on B
368  if (part == 0 && where == InKet) {
369  F a(Tint->bra(0,0));
370  F b(Tint->ket(0,0) - _1);
371  if (!exists(b)) return;
372  F c(Tint->bra(1,0));
373  F d(Tint->ket(1,0));
374 
375  auto ABCD_K = factory.make_child(a,b,c,d,0u,oK);
376  if (is_simple()) { expr_ = Vector("R12kG12_pfac0_0")[dir] * ABCD_K; nflops_+=1; }
377  const F& bm1 = b - _1;
378  if (exists(bm1)) {
379  auto ABm1CD_K = factory.make_child(a,bm1,c,d,0u,oK);
380  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("R12kG12_pfac1_0") * ABm1CD_K; nflops_+=3; }
381  }
382  const F& dm1 = d - _1;
383  if (exists(dm1)) {
384  auto ABCDm1_K = factory.make_child(a,b,c,dm1,0u,oK);
385  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("R12kG12_pfac2") * ABCDm1_K; nflops_+=3; }
386  }
387  if (K != 0) {
388  const R12kG12 oKm2((R12kG12Descr(K-2)));
389  auto ABp1CD_Km2 = factory.make_child(a,b+_1,c,d,0u,oKm2);
390  auto ABCDp1_Km2 = factory.make_child(a,b,c,d+_1,0u,oKm2);
391  auto ABCD_Km2 = factory.make_child(a,b,c,d,0u,oKm2);
392  if (is_simple()) { expr_ += Scalar((double)K) * Scalar("R12kG12_pfac3_0")
393  * (ABp1CD_Km2 - ABCDp1_Km2 + Vector("R12kG12_pfac4_0")[dir] * ABCD_Km2); nflops_+=6; }
394  }
395  return;
396  }
397  // Build on D
398  if (part == 1 && where == InKet) {
399  F a(Tint->bra(0,0));
400  F b(Tint->ket(0,0));
401  F c(Tint->bra(1,0));
402  F d(Tint->ket(1,0) - _1);
403  if (!exists(d)) return;
404 
405  auto ABCD_K = factory.make_child(a,b,c,d,0u,oK);
406  if (is_simple()) { expr_ = Vector("R12kG12_pfac0_1")[dir] * ABCD_K; nflops_+=1; }
407  const F& dm1 = d - _1;
408  if (exists(dm1)) {
409  auto ABCDm1_K = factory.make_child(a,b,c,dm1,0u,oK);
410  if (is_simple()) { expr_ += Scalar(d[dir]) * Scalar("R12kG12_pfac1_1") * ABCDm1_K; nflops_+=3; }
411  }
412  const F& bm1 = b - _1;
413  if (exists(bm1)) {
414  auto ABm1CD_K = factory.make_child(a,bm1,c,d,0u,oK);
415  if (is_simple()) { expr_ += Scalar(b[dir]) * Scalar("R12kG12_pfac2") * ABm1CD_K; nflops_+=3; }
416  }
417  if (K != 0) {
418  const R12kG12 oKm2((R12kG12Descr(K-2)));
419  auto ABCDp1_Km2 = factory.make_child(a,b,c,d+_1,0u,oKm2);
420  auto ABp1CD_Km2 = factory.make_child(a,b+_1,c,d,0u,oKm2);
421  auto ABCD_Km2 = factory.make_child(a,b,c,d,0u,oKm2);
422  if (is_simple()) { expr_ += Scalar((double)K) * Scalar("R12kG12_pfac3_1")
423  * (ABCDp1_Km2 - ABp1CD_Km2 + Vector("R12kG12_pfac4_1")[dir] * ABCD_Km2); nflops_+=6; }
424  }
425  return;
426  }
427  } // end of 0b0d case
428 
429  return;
430  } // K >= 0
431  }
432 
433  #if LIBINT_ENABLE_GENERIC_CODE
434  template <class F, int part, FunctionPosition where>
435  bool
436  VRR_11_R12kG12_11<F,part,where>::has_generic(const SafePtr<CompilationParameters>& cparams) const
437  {
438  F sh_a(target_->bra(0,0));
439  F sh_b(target_->ket(0,0));
440  F sh_c(target_->bra(1,0));
441  F sh_d(target_->ket(1,0));
442  const unsigned int max_opt_am = cparams->max_am_opt();
443  // generic code works for a0c0 and 0a0c classes where am(a) > 1 and am(c) > 1
444  // to generate optimized code for xxxx integral need to generate specialized code for up to (x+x)0(x+x)0 integrals
446  sh_b.zero() && sh_d.zero() &&
447  (sh_a.norm() > std::max(2*max_opt_am,1u) ||
448  sh_c.norm() > std::max(2*max_opt_am,1u)
449  ) &&
450  (sh_a.norm() > 1u && sh_c.norm() > 1u)
451  )
452  return true;
454  sh_a.zero() && sh_c.zero() &&
455  (sh_b.norm() > std::max(2*max_opt_am,1u) ||
456  sh_d.norm() > std::max(2*max_opt_am,1u)
457  ) &&
458  (sh_b.norm() > 1u && sh_d.norm() > 1u)
459  )
460  return true;
461  return false;
462  }
463 
464  template <class F, int part, FunctionPosition where>
465  std::string
467  {
468  F sh_a(target_->bra(0,0));
469  F sh_b(target_->ket(0,0));
470  F sh_c(target_->bra(1,0));
471  F sh_d(target_->ket(1,0));
472  const bool xsxs = sh_b.zero() && sh_d.zero();
473  const bool sxsx = sh_a.zero() && sh_c.zero();
474  const int K = target_->oper()->descr().K();
475  if (K == -1) {
476  if (xsxs)
477  return std::string("OSVRR_xs_xs.h");
478  if (sxsx)
479  return std::string("OSVRR_sx_sx.h");
480  }
481  else {
482  return std::string("VRR_r12kg12_xs_xs.h");
483  }
484  assert(false); // unreachable
485  }
486 
487  template <class F, int part, FunctionPosition where>
488  std::string
489  VRR_11_R12kG12_11<F,part,where>::generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const
490  {
491  const int K = target_->oper()->descr().K();
492  std::ostringstream oss;
493  F sh_a(target_->bra(0,0));
494  F sh_b(target_->ket(0,0));
495  F sh_c(target_->bra(1,0));
496  F sh_d(target_->ket(1,0));
497  const bool xsxs = sh_b.zero() && sh_d.zero();
498  const bool sxsx = sh_a.zero() && sh_c.zero();
499 
500  bool ahlrichs_simplification = false;
501  bool unit_s = false;
502  if (xsxs) {
503  unit_s = sh_b.is_unit();
504  }
505  else { // (sxsx)
506  unit_s = sh_a.is_unit();
507  }
508 
509  oss << "using namespace libint2;" << endl;
510  if (K == -1) {
511  if (xsxs) {
512  oss << "libint2::OSVRR_xs_xs<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",";
513  if (not ahlrichs_simplification) oss << (unit_s ? "true," : "false,");
514  oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
515  }
516  if (sxsx) {
517  oss << "libint2::OSVRR_sx_sx<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",";
518  if (not ahlrichs_simplification) oss << (unit_s ? "true," : "false,");
519  oss << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
520  }
521  }
522  else {
523  if (xsxs) {
524  oss << "libint2::VRR_r12kg12_xs_xs<" << part << "," << sh_a.norm() << "," << sh_c.norm() << ",";
525  oss << K << "," << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
526  }
527  if (sxsx) {
528  // NOTE that using same function, the only difference is in the RR prefactors
529  oss << "libint2::VRR_r12kg12_xs_xs<" << part << "," << sh_b.norm() << "," << sh_d.norm() << ",";
530  oss << K << "," << ((context->cparams()->max_vector_length() == 1) ? "false" : "true");
531  }
532  }
533  oss << ">::compute(inteval";
534 
535  // if K == -1 follow the same logic as for standard VRR
536  if (K == -1) {
537  oss << "," << args->symbol(0); // target
538  if (not ahlrichs_simplification && unit_s) // purely to avoid having a 4-term generic RR, reuse 5-term with dummy argument
539  oss << ",0"; // src0-> 0x0
540  const unsigned int nargs = args->n();
541  for(unsigned int a=1; a<nargs; a++) {
542  oss << "," << args->symbol(a);
543  }
544  }
545  else { // K != -1
546  const unsigned int nargs = args->n();
547  for(unsigned int a=0; a<nargs; a++) {
548  oss << "," << args->symbol(a);
549  }
550  }
551  // if K == 0 add dummy arguments so that the same generic function can be used for all K>=0 cases
552  if (K == 0) {
553  for(unsigned int a=0; a<3; ++a) {
554  oss << ",(LIBINT2_REALTYPE*)0";
555  }
556  }
557 
558  oss << ");";
559 
560  return oss.str();
561  }
562  #endif // #if !LIBINT_ENABLE_GENERIC_CODE
563 
564 };
565 
566 #endif
R12_k_G12 is a two-body operator of form r_{12}^k * exp(- * r_{12}), where k is an integer and is a ...
Definition: oper.h:298
const SafePtr< DGVertex > & make_child(const F &A, const F &B, const F &C, const F &D, const AuxIndexType &aux=AuxIndexType(), const OperType &oper=OperType())
make_child
Definition: generic_rr.h:157
TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition: bfset.h:879
static bool directional()
Default directionality.
Definition: vrr_11_r12kg12_11.h:50
VRR Recurrence Relation for 2-e integrals of the R12_k_G12 operators.
Definition: vrr_11_r12kg12_11.h:35
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
Definition: stdarray.h:18
Definition: prefactors.h:159
GenOper is a single operator described by descriptor Descr.
Definition: oper.h:152
Generic integral over a two-body operator with one bfs for each particle in bra and ket...
Definition: integral_11_11.h:32
RRImpl must inherit GenericRecurrenceRelation<RRImpl>
Definition: generic_rr.h:48
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
Definition: algebra.h:32
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:91
Set of basis functions.
Definition: bfset.h:41
bool zero() const
norm() == 0
Definition: bfset.h:153
Helps GenericRecurrenceRelation to work around the compiler problem with make_child.
Definition: generic_rr.h:148