LIBINT 2.9.0
vrr_11_r12kg12_11.h
1/*
2 * Copyright (C) 2004-2024 Edward F. Valeev
3 *
4 * This file is part of Libint compiler.
5 *
6 * Libint compiler 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 compiler 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 compiler. 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
33template <class BFSet, int part, FunctionPosition where>
34class VRR_11_R12kG12_11 : public GenericRecurrenceRelation<
35 VRR_11_R12kG12_11<BFSet, part, where>, BFSet,
36 GenIntegralSet_11_11<BFSet, R12kG12, mType> > {
37 public:
38 typedef VRR_11_R12kG12_11 ThisType;
39 typedef BFSet BasisFunctionType;
40 typedef GenIntegralSet_11_11<BFSet, R12kG12, mType> TargetType;
41 typedef GenericRecurrenceRelation<ThisType, BFSet, TargetType> ParentType;
42 friend class GenericRecurrenceRelation<ThisType, BFSet, TargetType>;
43 static const unsigned int max_nchildren = 8;
44
45 using ParentType::Instance;
46
48 static bool directional() { return ParentType::default_directional(); }
49
50 private:
51 using ParentType::is_simple;
52 using ParentType::target_;
53 using ParentType::RecurrenceRelation::expr_;
54 using ParentType::RecurrenceRelation::nflops_;
55
58 VRR_11_R12kG12_11(const std::shared_ptr<TargetType>&, unsigned int dir);
59
60 static std::string descr() { return "VRR"; }
65 std::string generate_label() const override {
66 typedef typename TargetType::AuxIndexType mType;
67 static std::shared_ptr<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,
71 target_->oper());
72 return os.str();
73 }
74
75#if LIBINT_ENABLE_GENERIC_CODE
77 bool has_generic(
78 const std::shared_ptr<CompilationParameters>& cparams) const override;
80 std::string generic_header() const override;
82 std::string generic_instance(
83 const std::shared_ptr<CodeContext>& context,
84 const std::shared_ptr<CodeSymbols>& args) const override;
85#endif
86};
87
88template <class F, int part, FunctionPosition where>
89VRR_11_R12kG12_11<F, part, where>::VRR_11_R12kG12_11(
90 const std::shared_ptr<TargetType>& Tint, unsigned int dir)
91 : ParentType(Tint, dir) {
92 using namespace libint2::algebra;
93 using namespace libint2::prefactor;
94 const int K = Tint->oper()->descr().K();
95 typedef R12_k_G12_Descr R12kG12Descr;
96 const R12kG12 oK((R12kG12Descr(K)));
97 const unsigned int m = Tint->aux()->elem(0);
98 const F _1 = unit<F>(dir);
99
100 // does not work for contracted functions
101 {
102 F a(Tint->bra(0, 0));
103 F b(Tint->ket(0, 0));
104 F c(Tint->bra(1, 0));
105 F d(Tint->ket(1, 0));
106 if (a.contracted() || b.contracted() || c.contracted() || 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() || !dB.zero() || !dC.zero() || !dD.zero();
117 assert(!deriv);
118
119 typedef TargetType ChildType;
120 ChildFactory<ThisType, ChildType> factory(this);
121
122 // if K is -1, the recurrence relation looks exactly as it would for ERI
123 // thus generate the same code, and remember to use appropriate prefactors
124 if (K == -1) {
125 // Build on A
126 if (part == 0 && where == InBra) {
127 F a(Tint->bra(0, 0) - _1);
128 if (!exists(a)) return;
129 F b(Tint->ket(0, 0));
130 F c(Tint->bra(1, 0));
131 F d(Tint->ket(1, 0));
132
133 auto ABCD_m = factory.make_child(a, b, c, d, m, oK);
134 auto ABCD_mp1 = factory.make_child(a, b, c, d, m + 1, oK);
135 if (is_simple()) {
136 expr_ = Vector("PA")[dir] * ABCD_m + Vector("WP")[dir] * ABCD_mp1;
137 nflops_ += 3;
138 }
139
140 const F& am1 = a - _1;
141 if (exists(am1)) {
142 auto Am1BCD_m = factory.make_child(am1, b, c, d, m, oK);
143 auto Am1BCD_mp1 = factory.make_child(am1, b, c, d, m + 1, oK);
144 if (is_simple()) {
145 expr_ += Scalar(a[dir]) * Scalar("oo2z") *
146 (Am1BCD_m - Scalar("roz") * Am1BCD_mp1);
147 nflops_ += 5;
148 }
149 }
150 const F& bm1 = b - _1;
151 if (exists(bm1)) {
152 auto ABm1CD_m = factory.make_child(a, bm1, c, d, m, oK);
153 auto ABm1CD_mp1 = factory.make_child(a, bm1, c, d, m + 1, oK);
154 if (is_simple()) {
155 expr_ += Scalar(b[dir]) * Scalar("oo2z") *
156 (ABm1CD_m - Scalar("roz") * ABm1CD_mp1);
157 nflops_ += 5;
158 }
159 }
160 const F& cm1 = c - _1;
161 if (exists(cm1)) {
162 auto ABCm1D_mp1 = factory.make_child(a, b, cm1, d, m + 1, oK);
163 if (is_simple()) {
164 expr_ += Scalar(c[dir]) * Scalar("oo2ze") * ABCm1D_mp1;
165 nflops_ += 3;
166 }
167 }
168 const F& dm1 = d - _1;
169 if (exists(dm1)) {
170 auto ABCDm1_mp1 = factory.make_child(a, b, c, dm1, m + 1, oK);
171 if (is_simple()) {
172 expr_ += Scalar(d[dir]) * Scalar("oo2ze") * ABCDm1_mp1;
173 nflops_ += 3;
174 }
175 }
176 return;
177 }
178 // Build on A
179 if (part == 0 && where == InKet) {
180 F a(Tint->bra(0, 0));
181 F b(Tint->ket(0, 0) - _1);
182 if (!exists(b)) return;
183 F c(Tint->bra(1, 0));
184 F d(Tint->ket(1, 0));
185
186 auto ABCD_m = factory.make_child(a, b, c, d, m, oK);
187 auto ABCD_mp1 = factory.make_child(a, b, c, d, m + 1, oK);
188 if (is_simple()) {
189 expr_ = Vector("PB")[dir] * ABCD_m + Vector("WP")[dir] * ABCD_mp1;
190 nflops_ += 3;
191 }
192
193 const F& am1 = a - _1;
194 if (exists(am1)) {
195 auto Am1BCD_m = factory.make_child(am1, b, c, d, m, oK);
196 auto Am1BCD_mp1 = factory.make_child(am1, b, c, d, m + 1, oK);
197 if (is_simple()) {
198 expr_ += Scalar(a[dir]) * Scalar("oo2z") *
199 (Am1BCD_m - Scalar("roz") * Am1BCD_mp1);
200 nflops_ += 5;
201 }
202 }
203 const F& bm1 = b - _1;
204 if (exists(bm1)) {
205 auto ABm1CD_m = factory.make_child(a, bm1, c, d, m, oK);
206 auto ABm1CD_mp1 = factory.make_child(a, bm1, c, d, m + 1, oK);
207 if (is_simple()) {
208 expr_ += Scalar(b[dir]) * Scalar("oo2z") *
209 (ABm1CD_m - Scalar("roz") * ABm1CD_mp1);
210 nflops_ += 5;
211 }
212 }
213 const F& cm1 = c - _1;
214 if (exists(cm1)) {
215 auto ABCm1D_mp1 = factory.make_child(a, b, cm1, d, m + 1, oK);
216 if (is_simple()) {
217 expr_ += Scalar(c[dir]) * Scalar("oo2ze") * ABCm1D_mp1;
218 nflops_ += 3;
219 }
220 }
221 const F& dm1 = d - _1;
222 if (exists(dm1)) {
223 auto ABCDm1_mp1 = factory.make_child(a, b, c, dm1, m + 1, oK);
224 if (is_simple()) {
225 expr_ += Scalar(d[dir]) * Scalar("oo2ze") * ABCDm1_mp1;
226 nflops_ += 3;
227 }
228 }
229 return;
230 }
231 // Build on C
232 if (part == 1 && where == InBra) {
233 F a(Tint->bra(0, 0));
234 F b(Tint->ket(0, 0));
235 F c(Tint->bra(1, 0) - _1);
236 if (!exists(c)) return;
237 F d(Tint->ket(1, 0));
238
239 auto ABCD_m = factory.make_child(a, b, c, d, m, oK);
240 auto ABCD_mp1 = factory.make_child(a, b, c, d, m + 1, oK);
241 if (is_simple()) {
242 expr_ = Vector("QC")[dir] * ABCD_m + Vector("WQ")[dir] * ABCD_mp1;
243 nflops_ += 3;
244 }
245
246 const F& cm1 = c - _1;
247 if (exists(cm1)) {
248 auto ABCm1D_m = factory.make_child(a, b, cm1, d, m, oK);
249 auto ABCm1D_mp1 = factory.make_child(a, b, cm1, d, m + 1, oK);
250 if (is_simple()) {
251 expr_ += Scalar(c[dir]) * Scalar("oo2e") *
252 (ABCm1D_m - Scalar("roe") * ABCm1D_mp1);
253 nflops_ += 5;
254 }
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()) {
261 expr_ += Scalar(d[dir]) * Scalar("oo2e") *
262 (ABCDm1_m - Scalar("roe") * ABCDm1_mp1);
263 nflops_ += 5;
264 }
265 }
266 const F& am1 = a - _1;
267 if (exists(am1)) {
268 auto Am1BCD_mp1 = factory.make_child(am1, b, c, d, m + 1, oK);
269 if (is_simple()) {
270 expr_ += Scalar(a[dir]) * Scalar("oo2ze") * Am1BCD_mp1;
271 nflops_ += 3;
272 }
273 }
274 const F& bm1 = b - _1;
275 if (exists(bm1)) {
276 auto ABm1CD_mp1 = factory.make_child(a, bm1, c, d, m + 1, oK);
277 if (is_simple()) {
278 expr_ += Scalar(b[dir]) * Scalar("oo2ze") * ABm1CD_mp1;
279 nflops_ += 3;
280 }
281 }
282 return;
283 }
284 // Build on D
285 if (part == 1 && where == InKet) {
286 F a(Tint->bra(0, 0));
287 F b(Tint->ket(0, 0));
288 F c(Tint->bra(1, 0));
289 F d(Tint->ket(1, 0) - _1);
290 if (!exists(d)) return;
291
292 auto ABCD_m = factory.make_child(a, b, c, d, m, oK);
293 auto ABCD_mp1 = factory.make_child(a, b, c, d, m + 1, oK);
294 if (is_simple()) {
295 expr_ = Vector("QD")[dir] * ABCD_m + Vector("WQ")[dir] * ABCD_mp1;
296 nflops_ += 3;
297 }
298
299 const F& cm1 = c - _1;
300 if (exists(cm1)) {
301 auto ABCm1D_m = factory.make_child(a, b, cm1, d, m, oK);
302 auto ABCm1D_mp1 = factory.make_child(a, b, cm1, d, m + 1, oK);
303 if (is_simple()) {
304 expr_ += Scalar(c[dir]) * Scalar("oo2e") *
305 (ABCm1D_m - Scalar("roe") * ABCm1D_mp1);
306 nflops_ += 5;
307 }
308 }
309 const F& dm1 = d - _1;
310 if (exists(dm1)) {
311 auto ABCDm1_m = factory.make_child(a, b, c, dm1, m, oK);
312 auto ABCDm1_mp1 = factory.make_child(a, b, c, dm1, m + 1, oK);
313 if (is_simple()) {
314 expr_ += Scalar(d[dir]) * Scalar("oo2e") *
315 (ABCDm1_m - Scalar("roe") * ABCDm1_mp1);
316 nflops_ += 5;
317 }
318 }
319 const F& am1 = a - _1;
320 if (exists(am1)) {
321 auto Am1BCD_mp1 = factory.make_child(am1, b, c, d, m + 1, oK);
322 if (is_simple()) {
323 expr_ += Scalar(a[dir]) * Scalar("oo2ze") * Am1BCD_mp1;
324 nflops_ += 3;
325 }
326 }
327 const F& bm1 = b - _1;
328 if (exists(bm1)) {
329 auto ABm1CD_mp1 = factory.make_child(a, bm1, c, d, m + 1, oK);
330 if (is_simple()) {
331 expr_ += Scalar(b[dir]) * Scalar("oo2ze") * ABm1CD_mp1;
332 nflops_ += 3;
333 }
334 }
335 return;
336 }
337 return;
338 } // K == -1?
339 else {
340 // K != -1, the auxiliary quantum number is not used
341 if (m != 0)
342 throw std::logic_error(
343 "VRR_11_R12kG12_11<I,F,K,part,where>::children_and_expr_Kge0() -- "
344 "nonzero auxiliary quantum detected.");
345
346 // can build (a0|c0) or (0b|0d)
347 bool xsxs = false;
348 bool sxsx = false;
349 {
350 F b(Tint->ket(0, 0) - _1);
351 F d(Tint->ket(1, 0) - _1);
352 if (!exists(b) && !exists(d)) xsxs = true;
353 }
354 {
355 F a(Tint->bra(0, 0) - _1);
356 F c(Tint->bra(1, 0) - _1);
357 if (!exists(a) && !exists(c)) sxsx = true;
358 }
359 // can't handle the general case
360 if (!xsxs && !sxsx) return;
361 // can't handle (ss|ss) case
362 if (xsxs && sxsx) return;
363
364 if (xsxs) {
365 // Build on A
366 if (part == 0 && where == InBra) {
367 F a(Tint->bra(0, 0) - _1);
368 if (!exists(a)) return;
369 F b(Tint->ket(0, 0));
370 F c(Tint->bra(1, 0));
371 F d(Tint->ket(1, 0));
372
373 auto ABCD_K = factory.make_child(a, b, c, d, 0u, oK);
374 if (is_simple()) {
375 expr_ = Vector("R12kG12_pfac0_0")[dir] * ABCD_K;
376 nflops_ += 1;
377 }
378 const F& am1 = a - _1;
379 if (exists(am1)) {
380 auto Am1BCD_K = factory.make_child(am1, b, c, d, 0u, oK);
381 if (is_simple()) {
382 expr_ += Scalar(a[dir]) * Scalar("R12kG12_pfac1_0") * Am1BCD_K;
383 nflops_ += 3;
384 }
385 }
386 const F& cm1 = c - _1;
387 if (exists(cm1)) {
388 auto ABCm1D_K = factory.make_child(a, b, cm1, d, 0u, oK);
389 if (is_simple()) {
390 expr_ += Scalar(c[dir]) * Scalar("R12kG12_pfac2") * ABCm1D_K;
391 nflops_ += 3;
392 }
393 }
394 if (K != 0) {
395 const R12kG12 oKm2((R12kG12Descr(K - 2)));
396 auto Ap1BCD_Km2 = factory.make_child(a + _1, b, c, d, 0u, oKm2);
397 auto ABCp1D_Km2 = factory.make_child(a, b, c + _1, d, 0u, oKm2);
398 auto ABCD_Km2 = factory.make_child(a, b, c, d, 0u, oKm2);
399 if (is_simple()) {
400 expr_ += Scalar((double)K) * Scalar("R12kG12_pfac3_0") *
401 (Ap1BCD_Km2 - ABCp1D_Km2 +
402 Vector("R12kG12_pfac4_0")[dir] * ABCD_Km2);
403 nflops_ += 6;
404 }
405 }
406 return;
407 }
408 // Build on C
409 if (part == 1 && where == InBra) {
410 F a(Tint->bra(0, 0));
411 F b(Tint->ket(0, 0));
412 F c(Tint->bra(1, 0) - _1);
413 if (!exists(c)) return;
414 F d(Tint->ket(1, 0));
415
416 auto ABCD_K = factory.make_child(a, b, c, d, 0u, oK);
417 if (is_simple()) {
418 expr_ = Vector("R12kG12_pfac0_1")[dir] * ABCD_K;
419 nflops_ += 1;
420 }
421 const F& cm1 = c - _1;
422 if (exists(cm1)) {
423 auto ABCm1D_K = factory.make_child(a, b, cm1, d, 0u, oK);
424 if (is_simple()) {
425 expr_ += Scalar(c[dir]) * Scalar("R12kG12_pfac1_1") * ABCm1D_K;
426 nflops_ += 3;
427 }
428 }
429 const F& am1 = a - _1;
430 if (exists(am1)) {
431 auto Am1BCD_K = factory.make_child(am1, b, c, d, 0u, oK);
432 if (is_simple()) {
433 expr_ += Scalar(a[dir]) * Scalar("R12kG12_pfac2") * Am1BCD_K;
434 nflops_ += 3;
435 }
436 }
437 if (K != 0) {
438 const R12kG12 oKm2((R12kG12Descr(K - 2)));
439 auto ABCp1D_Km2 = factory.make_child(a, b, c + _1, d, 0u, oKm2);
440 auto Ap1BCD_Km2 = factory.make_child(a + _1, b, c, d, 0u, oKm2);
441 auto ABCD_Km2 = factory.make_child(a, b, c, d, 0u, oKm2);
442 if (is_simple()) {
443 expr_ += Scalar((double)K) * Scalar("R12kG12_pfac3_1") *
444 (ABCp1D_Km2 - Ap1BCD_Km2 +
445 Vector("R12kG12_pfac4_1")[dir] * ABCD_Km2);
446 nflops_ += 6;
447 }
448 }
449 return;
450 }
451 } // end of a0c0 case
452
453 if (sxsx) {
454 // Build on B
455 if (part == 0 && where == InKet) {
456 F a(Tint->bra(0, 0));
457 F b(Tint->ket(0, 0) - _1);
458 if (!exists(b)) return;
459 F c(Tint->bra(1, 0));
460 F d(Tint->ket(1, 0));
461
462 auto ABCD_K = factory.make_child(a, b, c, d, 0u, oK);
463 if (is_simple()) {
464 expr_ = Vector("R12kG12_pfac0_0")[dir] * ABCD_K;
465 nflops_ += 1;
466 }
467 const F& bm1 = b - _1;
468 if (exists(bm1)) {
469 auto ABm1CD_K = factory.make_child(a, bm1, c, d, 0u, oK);
470 if (is_simple()) {
471 expr_ += Scalar(b[dir]) * Scalar("R12kG12_pfac1_0") * ABm1CD_K;
472 nflops_ += 3;
473 }
474 }
475 const F& dm1 = d - _1;
476 if (exists(dm1)) {
477 auto ABCDm1_K = factory.make_child(a, b, c, dm1, 0u, oK);
478 if (is_simple()) {
479 expr_ += Scalar(d[dir]) * Scalar("R12kG12_pfac2") * ABCDm1_K;
480 nflops_ += 3;
481 }
482 }
483 if (K != 0) {
484 const R12kG12 oKm2((R12kG12Descr(K - 2)));
485 auto ABp1CD_Km2 = factory.make_child(a, b + _1, c, d, 0u, oKm2);
486 auto ABCDp1_Km2 = factory.make_child(a, b, c, d + _1, 0u, oKm2);
487 auto ABCD_Km2 = factory.make_child(a, b, c, d, 0u, oKm2);
488 if (is_simple()) {
489 expr_ += Scalar((double)K) * Scalar("R12kG12_pfac3_0") *
490 (ABp1CD_Km2 - ABCDp1_Km2 +
491 Vector("R12kG12_pfac4_0")[dir] * ABCD_Km2);
492 nflops_ += 6;
493 }
494 }
495 return;
496 }
497 // Build on D
498 if (part == 1 && where == InKet) {
499 F a(Tint->bra(0, 0));
500 F b(Tint->ket(0, 0));
501 F c(Tint->bra(1, 0));
502 F d(Tint->ket(1, 0) - _1);
503 if (!exists(d)) return;
504
505 auto ABCD_K = factory.make_child(a, b, c, d, 0u, oK);
506 if (is_simple()) {
507 expr_ = Vector("R12kG12_pfac0_1")[dir] * ABCD_K;
508 nflops_ += 1;
509 }
510 const F& dm1 = d - _1;
511 if (exists(dm1)) {
512 auto ABCDm1_K = factory.make_child(a, b, c, dm1, 0u, oK);
513 if (is_simple()) {
514 expr_ += Scalar(d[dir]) * Scalar("R12kG12_pfac1_1") * ABCDm1_K;
515 nflops_ += 3;
516 }
517 }
518 const F& bm1 = b - _1;
519 if (exists(bm1)) {
520 auto ABm1CD_K = factory.make_child(a, bm1, c, d, 0u, oK);
521 if (is_simple()) {
522 expr_ += Scalar(b[dir]) * Scalar("R12kG12_pfac2") * ABm1CD_K;
523 nflops_ += 3;
524 }
525 }
526 if (K != 0) {
527 const R12kG12 oKm2((R12kG12Descr(K - 2)));
528 auto ABCDp1_Km2 = factory.make_child(a, b, c, d + _1, 0u, oKm2);
529 auto ABp1CD_Km2 = factory.make_child(a, b + _1, c, d, 0u, oKm2);
530 auto ABCD_Km2 = factory.make_child(a, b, c, d, 0u, oKm2);
531 if (is_simple()) {
532 expr_ += Scalar((double)K) * Scalar("R12kG12_pfac3_1") *
533 (ABCDp1_Km2 - ABp1CD_Km2 +
534 Vector("R12kG12_pfac4_1")[dir] * ABCD_Km2);
535 nflops_ += 6;
536 }
537 }
538 return;
539 }
540 } // end of 0b0d case
541
542 return;
543 } // K >= 0
544}
545
546#if LIBINT_ENABLE_GENERIC_CODE
547template <class F, int part, FunctionPosition where>
548bool VRR_11_R12kG12_11<F, part, where>::has_generic(
549 const std::shared_ptr<CompilationParameters>& cparams) const {
550 F sh_a(target_->bra(0, 0));
551 F sh_b(target_->ket(0, 0));
552 F sh_c(target_->bra(1, 0));
553 F sh_d(target_->ket(1, 0));
554 const unsigned int max_opt_am = cparams->max_am_opt();
555 // generic code works for a0c0 and 0a0c classes where am(a) > 1 and am(c) > 1
556 // to generate optimized code for xxxx integral need to generate specialized
557 // code for up to (x+x)0(x+x)0 integrals
558 if (!TrivialBFSet<F>::result && sh_b.zero() && sh_d.zero() &&
559 (sh_a.norm() > std::max(2 * max_opt_am, 1u) ||
560 sh_c.norm() > std::max(2 * max_opt_am, 1u)) &&
561 (sh_a.norm() > 1u && sh_c.norm() > 1u))
562 return true;
563 if (!TrivialBFSet<F>::result && sh_a.zero() && sh_c.zero() &&
564 (sh_b.norm() > std::max(2 * max_opt_am, 1u) ||
565 sh_d.norm() > std::max(2 * max_opt_am, 1u)) &&
566 (sh_b.norm() > 1u && sh_d.norm() > 1u))
567 return true;
568 return false;
569}
570
571template <class F, int part, FunctionPosition where>
572std::string VRR_11_R12kG12_11<F, part, where>::generic_header() const {
573 F sh_a(target_->bra(0, 0));
574 F sh_b(target_->ket(0, 0));
575 F sh_c(target_->bra(1, 0));
576 F sh_d(target_->ket(1, 0));
577 const bool xsxs = sh_b.zero() && sh_d.zero();
578 const bool sxsx = sh_a.zero() && sh_c.zero();
579 const int K = target_->oper()->descr().K();
580 if (K == -1) {
581 if (xsxs) return std::string("OSVRR_xs_xs.h");
582 if (sxsx) return std::string("OSVRR_sx_sx.h");
583 } else {
584 return std::string("VRR_r12kg12_xs_xs.h");
585 }
586 assert(false); // unreachable
587}
588
589template <class F, int part, FunctionPosition where>
590std::string VRR_11_R12kG12_11<F, part, where>::generic_instance(
591 const std::shared_ptr<CodeContext>& context,
592 const std::shared_ptr<CodeSymbols>& args) const {
593 const int K = target_->oper()->descr().K();
594 std::ostringstream oss;
595 F sh_a(target_->bra(0, 0));
596 F sh_b(target_->ket(0, 0));
597 F sh_c(target_->bra(1, 0));
598 F sh_d(target_->ket(1, 0));
599 const bool xsxs = sh_b.zero() && sh_d.zero();
600 const bool sxsx = sh_a.zero() && sh_c.zero();
601
602 bool ahlrichs_simplification = false;
603 bool unit_s = false;
604 if (xsxs) {
605 unit_s = sh_b.is_unit();
606 } else { // (sxsx)
607 unit_s = sh_a.is_unit();
608 }
609
610 oss << "using namespace libint2;" << std::endl;
611 if (K == -1) {
612 if (xsxs) {
613 oss << "libint2::OSVRR_xs_xs<" << part << "," << sh_a.norm() << ","
614 << sh_c.norm() << ",";
615 if (not ahlrichs_simplification) oss << (unit_s ? "true," : "false,");
616 oss << ((context->cparams()->max_vector_length() == 1) ? "false"
617 : "true");
618 }
619 if (sxsx) {
620 oss << "libint2::OSVRR_sx_sx<" << part << "," << sh_b.norm() << ","
621 << sh_d.norm() << ",";
622 if (not ahlrichs_simplification) oss << (unit_s ? "true," : "false,");
623 oss << ((context->cparams()->max_vector_length() == 1) ? "false"
624 : "true");
625 }
626 } else {
627 if (xsxs) {
628 oss << "libint2::VRR_r12kg12_xs_xs<" << part << "," << sh_a.norm() << ","
629 << sh_c.norm() << ",";
630 oss << K << ","
631 << ((context->cparams()->max_vector_length() == 1) ? "false"
632 : "true");
633 }
634 if (sxsx) {
635 // NOTE that using same function, the only difference is in the RR
636 // prefactors
637 oss << "libint2::VRR_r12kg12_xs_xs<" << part << "," << sh_b.norm() << ","
638 << sh_d.norm() << ",";
639 oss << K << ","
640 << ((context->cparams()->max_vector_length() == 1) ? "false"
641 : "true");
642 }
643 }
644 oss << ">::compute(inteval";
645
646 // if K == -1 follow the same logic as for standard VRR
647 if (K == -1) {
648 oss << "," << args->symbol(0); // target
649 if (not ahlrichs_simplification &&
650 unit_s) // purely to avoid having a 4-term generic RR, reuse 5-term
651 // with dummy argument
652 oss << ",0"; // src0-> 0x0
653 const unsigned int nargs = args->n();
654 for (unsigned int a = 1; a < nargs; a++) {
655 oss << "," << args->symbol(a);
656 }
657 } else { // K != -1
658 const unsigned int nargs = args->n();
659 for (unsigned int a = 0; a < nargs; a++) {
660 oss << "," << args->symbol(a);
661 }
662 }
663 // if K == 0 add dummy arguments so that the same generic function can be used
664 // for all K>=0 cases
665 if (K == 0) {
666 for (unsigned int a = 0; a < 3; ++a) {
667 oss << ",(LIBINT2_REALTYPE*)0";
668 }
669 }
670
671 oss << ");";
672
673 return oss.str();
674}
675#endif // #if !LIBINT_ENABLE_GENERIC_CODE
676
677}; // namespace libint2
678
679#endif
libint2::BFSet
Set of basis functions.
Definition bfset.h:43
libint2::ChildFactory
Helps GenericRecurrenceRelation to work around the compiler problem with make_child.
Definition generic_rr.h:150
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:36
libint2::GenOper
GenOper is a single operator described by descriptor Descr.
Definition oper.h:164
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:81
libint2::GenericRecurrenceRelation::default_directional
static bool default_directional()
is this recurrence relation parameterized by a direction (x, y, or z).
Definition generic_rr.h:101
libint2::GenericRecurrenceRelation::Instance
static std::shared_ptr< RRImpl > Instance(const std::shared_ptr< TargetType > &Tint, unsigned int dir)
Return an instance if applicable, or a null pointer otherwise.
Definition generic_rr.h:55
libint2::OriginDerivative
Represents cartesian derivatives of atom-centered basis functions.
Definition bfset.h:101
libint2::OriginDerivative::zero
bool zero() const
norm() == 0
Definition bfset.h:158
libint2::R12_k_G12_Descr
R12_k_G12 is a two-body operator of form r_{12}^k * exp(-\gamma * r_{12}), where k is an integer and ...
Definition oper.h:418
libint2::VRR_11_R12kG12_11
VRR Recurrence Relation for 2-e integrals of the R12_k_G12 operators.
Definition vrr_11_r12kg12_11.h:36
libint2::VRR_11_R12kG12_11::directional
static bool directional()
Default directionality.
Definition vrr_11_r12kg12_11.h:48
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:96
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:395
libint2::to_string
std::string to_string(const T &x)
Converts x to its string representation.
Definition entity.h:121
libint2::TrivialBFSet
TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition bfset.h:906
Generated by doxygen 1.10.0