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