LIBINT 2.9.0
engine.impl.h
1/*
2 * Copyright (C) 2004-2024 Edward F. Valeev
3 *
4 * This file is part of Libint library.
5 *
6 * Libint library is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser 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 library 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 Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with Libint library. If not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21#ifndef _libint2_src_lib_libint_engineimpl_h_
22#define _libint2_src_lib_libint_engineimpl_h_
23
24#include <iterator>
25
26#include "./deriv_map.h"
27#include "./engine.h"
28
29#pragma GCC diagnostic push
30#pragma GCC system_header
31#include <Eigen/Core>
32#pragma GCC diagnostic pop
33
34#include <libint2/boys.h>
35#if LIBINT_HAS_SYSTEM_BOOST_PREPROCESSOR_VARIADICS
36#include <boost/preprocessor.hpp>
37#include <boost/preprocessor/facilities/is_1.hpp>
38#else // use bundled boost
39#include <libint2/boost/preprocessor.hpp>
40#include <libint2/boost/preprocessor/facilities/is_1.hpp>
41#endif
42
43// extra PP macros
44
45#define BOOST_PP_MAKE_TUPLE_INTERNAL(z, i, last) \
46 i BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, last))
48#define BOOST_PP_MAKE_TUPLE(n) \
49 (BOOST_PP_REPEAT(n, BOOST_PP_MAKE_TUPLE_INTERNAL, BOOST_PP_DEC(n)))
50
51// the engine will be profiled by default if library was configured with
52// --enable-profile
53#ifdef LIBINT2_PROFILE
54#define LIBINT2_ENGINE_TIMERS
55// uncomment if want to profile each integral class
56#define LIBINT2_ENGINE_PROFILE_CLASS
57#endif
58// uncomment if want to profile the engine even if library was configured
59// without --enable-profile
60// # define LIBINT2_ENGINE_TIMERS
61
62namespace libint2 {
63
64template <typename T, unsigned N>
65typename std::remove_all_extents<T>::type* to_ptr1(T (&a)[N]) {
66 return reinterpret_cast<typename std::remove_all_extents<T>::type*>(&a);
67}
68
73#define BOOST_PP_NBODY_OPERATOR_LIST \
74 (overlap, \
75 (kinetic, \
76 (elecpot, \
77 (elecpot, \
78 (elecpot, \
79 (1emultipole, \
80 (2emultipole, \
81 (3emultipole, \
82 (sphemultipole, \
83 (opVop, \
84 (eri, \
85 (eri, \
86 (eri, \
87 (eri, \
88 (eri, \
89 (eri, \
90 (eri, (eri, (eri, (eri, BOOST_PP_NIL))))))))))))))))))))
91
92#define BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE \
93 BOOST_PP_MAKE_TUPLE(BOOST_PP_LIST_SIZE(BOOST_PP_NBODY_OPERATOR_LIST))
94#define BOOST_PP_NBODY_OPERATOR_INDEX_LIST \
95 BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE)
96#define BOOST_PP_NBODY_OPERATOR_LAST_ONEBODY_INDEX \
97 9 // opVop, the 10th member of BOOST_PP_NBODY_OPERATOR_LIST, is the last
98 // 1-body operator
99
100// make list of braket indices for n-body ints
101#define BOOST_PP_NBODY_BRAKET_INDEX_TUPLE \
102 BOOST_PP_MAKE_TUPLE(BOOST_PP_INC(BOOST_PP_NBODY_BRAKET_MAX_INDEX))
103#define BOOST_PP_NBODY_BRAKET_INDEX_LIST \
104 BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_BRAKET_INDEX_TUPLE)
105#define BOOST_PP_NBODY_BRAKET_RANK_TUPLE (2, 3, 4)
106#define BOOST_PP_NBODY_BRAKET_RANK_LIST \
107 BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_BRAKET_RANK_TUPLE)
108
109// make list of derivative orders for n-body ints
110#define BOOST_PP_NBODY_DERIV_ORDER_TUPLE \
111 BOOST_PP_MAKE_TUPLE(BOOST_PP_INC(LIBINT2_MAX_DERIV_ORDER))
112#define BOOST_PP_NBODY_DERIV_ORDER_LIST \
113 BOOST_PP_TUPLE_TO_LIST(BOOST_PP_NBODY_DERIV_ORDER_TUPLE)
114
116__libint2_engine_inline libint2::any default_params(const Operator& oper) {
117 switch (static_cast<int>(oper)) {
118#define BOOST_PP_NBODYENGINE_MCR1(r, data, i, elem) \
119 case i: \
120 return operator_traits<static_cast<Operator>(i)>::default_params();
121 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR1, _,
122 BOOST_PP_NBODY_OPERATOR_LIST)
123 default:
124 break;
125 }
126 assert(false && "missing case in switch"); // unreachable
127 abort();
128}
129
131
139template <typename... ShellPack>
140__libint2_engine_inline const Engine::target_ptr_vec& Engine::compute(
141 const libint2::Shell& first_shell, const ShellPack&... rest_of_shells) {
142 constexpr auto nargs = 1 + sizeof...(rest_of_shells);
143 assert(nargs == braket_rank() &&
144 "# of arguments to compute() does not match the braket type");
145
146 std::array<std::reference_wrapper<const Shell>, nargs> shells{
147 {first_shell, rest_of_shells...}};
148
149 if (operator_rank() == 1) {
150 if (nargs == 2) return compute1(shells[0], shells[1]);
151 } else if (operator_rank() == 2) {
152 auto compute_ptr_idx = ((static_cast<int>(oper_) -
153 static_cast<int>(Operator::first_2body_oper)) *
154 nbrakets_2body +
155 (static_cast<int>(braket_) -
156 static_cast<int>(BraKet::first_2body_braket))) *
157 nderivorders_2body +
158 deriv_order_;
159 assert(compute_ptr_idx >= 0 && compute_ptr_idx < compute2_ptrs().size());
160 auto compute_ptr = compute2_ptrs()[compute_ptr_idx];
161 assert(compute_ptr != nullptr && "2-body compute function not found");
162 if (nargs == 2)
163 return (this->*compute_ptr)(shells[0], Shell::unit(), shells[1],
164 Shell::unit(), nullptr, nullptr);
165 if (nargs == 3)
166 return (this->*compute_ptr)(shells[0], Shell::unit(), shells[1],
167 shells[2], nullptr, nullptr);
168 if (nargs == 4)
169 return (this->*compute_ptr)(shells[0], shells[1], shells[2], shells[3],
170 nullptr, nullptr);
171 }
172
173 assert(false && "missing feature"); // only reached if missing a feature
174 abort();
175}
176
181__libint2_engine_inline const Engine::target_ptr_vec& Engine::compute1(
182 const libint2::Shell& s1, const libint2::Shell& s2) {
183 // can only handle 1 contraction at a time
184 assert((s1.ncontr() == 1 && s2.ncontr() == 1) &&
185 "generally-contracted shells not yet supported");
186
187 const auto oper_is_nuclear =
188 (oper_ == Operator::nuclear || oper_ == Operator::erf_nuclear ||
189 oper_ == Operator::erfc_nuclear || oper_ == Operator::opVop);
190
191 const auto l1 = s1.contr[0].l;
192 const auto l2 = s2.contr[0].l;
193 assert(l1 <= lmax_ && "the angular momentum limit is exceeded");
194 assert(l2 <= lmax_ && "the angular momentum limit is exceeded");
195
196 // if want nuclear, make sure there is at least one nucleus .. otherwise the
197 // user likely forgot to call set_params
198 if (oper_is_nuclear && nparams() == 0)
199 throw std::logic_error(
200 "Engine requires charges, but no charges found; forgot to call "
201 "set_params()?");
202
203 const auto n1 = s1.size();
204 const auto n2 = s2.size();
205 const auto n12 = n1 * n2;
206 const auto ncart1 = s1.cartesian_size();
207 const auto ncart2 = s2.cartesian_size();
208 const auto ncart12 = ncart1 * ncart2;
209
210 // assert # of primitive pairs
211 const auto nprim1 = s1.nprim();
212 const auto nprim2 = s2.nprim();
213 const auto nprimpairs = nprim1 * nprim2;
214 assert(nprimpairs <= primdata_.size() &&
215 "the max number of primitive pairs exceeded");
216
217 auto nparam_sets = nparams();
218
219 // keep track if need to set targets_ explicitly
220 bool set_targets = set_targets_;
221
222 // # of targets computed by libint
223 const auto ntargets = nopers() * num_geometrical_derivatives(2, deriv_order_);
224
225 // Libint computes derivatives with respect to basis functions only, must
226 // must use translational invariance to recover derivatives w.r.t. operator
227 // degrees of freedom
228 // will compute derivs w.r.t. 2 Gaussian centers + (if nuclear) nparam_sets
229 // operator centers
230 const auto nderivcenters_shset = 2 + (oper_is_nuclear ? nparam_sets : 0);
231 const auto nderivcoord = 3 * nderivcenters_shset;
232 const auto num_shellsets_computed =
233 nopers() * num_geometrical_derivatives(nderivcenters_shset, deriv_order_);
234
235 // will use scratch_ if:
236 // - Coulomb ints are computed 1 charge at a time, contributions are
237 // accumulated in scratch_ (unless la==lb==0)
238 // - derivatives on the missing center need to be reconstructed (no need to
239 // accumulate into scratch though)
240 // NB ints in scratch are packed in order
241 const auto accumulate_ints_in_scratch = oper_is_nuclear;
242
243 // adjust max angular momentum, if needed
244 const auto lmax = std::max(l1, l2);
245 assert(lmax <= lmax_ && "the angular momentum limit is exceeded");
246
247 // N.B. for l=0 no need to transform to solid harmonics
248 // this is a workaround for the corner case of oper_ == Operator::*nuclear,
249 // and solid harmonics (s|s) integral ... beware the integral storage state
250 // machine
251 const auto tform_to_solids =
252 (s1.contr[0].pure || s2.contr[0].pure) && lmax != 0;
253
254 // simple (s|s) ints will be computed directly and accumulated in the first
255 // element of stack
256 const auto compute_directly =
257 lmax == 0 && deriv_order_ == 0 &&
258 (oper_ == Operator::overlap || oper_is_nuclear) &&
259 oper_ != Operator::opVop;
260 if (compute_directly) {
261 primdata_[0].stack[0] = 0;
262 targets_[0] = primdata_[0].stack;
263 }
264
265 if (accumulate_ints_in_scratch)
266 std::fill(std::begin(scratch_),
267 std::begin(scratch_) + num_shellsets_computed * ncart12, 0.0);
268
269 // loop over accumulation batches
270 for (auto pset = 0u; pset != nparam_sets; ++pset) {
271 if (!oper_is_nuclear)
272 assert(nparam_sets == 1 && "unexpected number of operator parameters");
273
274 auto p12 = 0;
275 for (auto p1 = 0; p1 != nprim1; ++p1) {
276 for (auto p2 = 0; p2 != nprim2; ++p2, ++p12) {
277 compute_primdata(primdata_[p12], s1, s2, p1, p2, pset);
278 }
279 }
280 primdata_[0].contrdepth = p12;
281
282 if (compute_directly) {
283 auto& result = primdata_[0].stack[0];
284 switch (oper_) {
285 case Operator::overlap:
286 for (auto p12 = 0; p12 != primdata_[0].contrdepth; ++p12)
287 result += primdata_[p12]._0_Overlap_0_x[0] *
288 primdata_[p12]._0_Overlap_0_y[0] *
289 primdata_[p12]._0_Overlap_0_z[0];
290 break;
291 case Operator::nuclear:
292 case Operator::erf_nuclear:
293 case Operator::erfc_nuclear:
294 for (auto p12 = 0; p12 != primdata_[0].contrdepth; ++p12)
295 result += primdata_[p12].LIBINT_T_S_ELECPOT_S(0)[0];
296 break;
297 default:
298 assert(false && "missing case in switch");
299 }
300 primdata_[0].targets[0] = &result;
301 } else {
302 const auto buildfnidx = s1.contr[0].l * hard_lmax_ + s2.contr[0].l;
303 assert(buildfnptrs_[buildfnidx] && "null build function ptr");
304 buildfnptrs_[buildfnidx](&primdata_[0]);
305
306 if (accumulate_ints_in_scratch) {
307 set_targets = true;
308 // - for non-derivative ints and first derivative ints the target
309 // ints computed by libint will appear at the front of targets_
310 // - for second and higher derivs need to re-index targets, hence
311 // will accumulate later, when computing operator derivatives via
312 // transinv
313 if (deriv_order_ <= 1) {
314 // accumulate targets computed by libint for this pset into the
315 // accumulated targets in scratch
316 auto s_target = &scratch_[0];
317 for (auto s = 0; s != ntargets; ++s, s_target += ncart12)
318 if (pset != 0)
319 std::transform(primdata_[0].targets[s],
320 primdata_[0].targets[s] + ncart12, s_target,
321 s_target, std::plus<value_type>());
322 else
323 std::copy(primdata_[0].targets[s],
324 primdata_[0].targets[s] + ncart12, s_target);
325 }
326
327 // 2. reconstruct derivatives of nuclear ints for each nucleus
328 // using translational invariance
329 // NB this is done in cartesian basis, otherwise would have to tform
330 // to solids contributions from every atom, rather than the running
331 // total at the end
332 if (deriv_order_ > 0) {
333 switch (deriv_order_) {
334 case 1: {
335 // first 6 shellsets are derivatives with respect to Gaussian
336 // positions
337 // following them are derivs with respect to nuclear coordinates
338 // (3 per nucleus)
339 assert(ntargets == 6 && "unexpected # of targets");
340 auto dest = &scratch_[0] + (6 + pset * 3) * ncart12;
341 for (auto s = 0; s != 3; ++s, dest += ncart12) {
342 auto src = primdata_[0].targets[s];
343 for (auto i = 0; i != ncart12; ++i) {
344 dest[i] = -src[i];
345 }
346 }
347 dest -= 3 * ncart12;
348 for (auto s = 3; s != 6; ++s, dest += ncart12) {
349 auto src = primdata_[0].targets[s];
350 for (auto i = 0; i != ncart12; ++i) {
351 dest[i] -= src[i];
352 }
353 }
354 } break;
355
356 case 2: {
357 // computes upper triangle index
358 // n2 = matrix size times 2
359 // i,j = indices, i<j
360 auto upper_triangle_index_ord = [](int n2, int i, int j) {
361 return i * (n2 - i - 1) / 2 + j;
362 };
363 // same as above, but orders i and j
364 auto upper_triangle_index = [&](int n2, int i, int j) {
365 return upper_triangle_index_ord(n2, std::min(i, j),
366 std::max(i, j));
367 };
368
369 // accumulate ints for this pset to scratch in locations
370 // remapped to overall deriv index
371 const auto ncoords_times_two = nderivcoord * 2;
372 for (auto d0 = 0, d01 = 0; d0 != 6; ++d0) {
373 for (auto d1 = d0; d1 != 6; ++d1, ++d01) {
374 const auto d01_full =
375 upper_triangle_index_ord(ncoords_times_two, d0, d1);
376 auto tgt = &scratch_[d01_full * ncart12];
377 if (pset != 0)
378 std::transform(primdata_[0].targets[d01],
379 primdata_[0].targets[d01] + ncart12, tgt,
380 tgt, std::plus<value_type>());
381 else
382 std::copy(primdata_[0].targets[d01],
383 primdata_[0].targets[d01] + ncart12, tgt);
384 }
385 }
386
387 // use translational invariance to build derivatives w.r.t.
388 // operator centers
389 {
390 // mixed derivatives: first deriv w.r.t. Gaussian, second
391 // w.r.t. operator coord pset
392 const auto c1 = 2 + pset;
393 for (auto c0 = 0; c0 != 2; ++c0) {
394 for (auto xyz0 = 0; xyz0 != 3; ++xyz0) {
395 const auto coord0 = c0 * 3 + xyz0;
396 for (auto xyz1 = 0; xyz1 != 3; ++xyz1) {
397 const auto coord1 = c1 * 3 + xyz1; // coord1 > coord0
398
399 const auto coord01_abs = upper_triangle_index_ord(
400 ncoords_times_two, coord0, coord1);
401 auto tgt = &scratch_[coord01_abs * ncart12];
402
403 // d2 / dAi dOj = - d2 / dAi dAj
404 {
405 auto coord1_A = xyz1;
406 const auto coord01_A =
407 upper_triangle_index(12, coord0, coord1_A);
408 const auto src = primdata_[0].targets[coord01_A];
409 for (auto i = 0; i != ncart12; ++i) tgt[i] = -src[i];
410 }
411
412 // d2 / dAi dOj -= d2 / dAi dBj
413 {
414 auto coord1_B = 3 + xyz1;
415 const auto coord01_B =
416 upper_triangle_index(12, coord0, coord1_B);
417 const auto src = primdata_[0].targets[coord01_B];
418 for (auto i = 0; i != ncart12; ++i) tgt[i] -= src[i];
419 }
420 }
421 }
422 }
423 } // mixed derivs
424 {
425 // operator derivs
426 const auto c0 = 2 + pset;
427 const auto c1 = c0;
428 for (auto xyz0 = 0; xyz0 != 3; ++xyz0) {
429 const auto coord0 = c0 * 3 + xyz0;
430 for (auto xyz1 = xyz0; xyz1 != 3; ++xyz1) {
431 const auto coord1 = c1 * 3 + xyz1; // coord1 > coord0
432
433 const auto coord01_abs = upper_triangle_index_ord(
434 ncoords_times_two, coord0, coord1);
435 auto tgt = &scratch_[coord01_abs * ncart12];
436
437 // d2 / dOi dOj = d2 / dAi dAj
438 {
439 auto coord0_A = xyz0;
440 auto coord1_A = xyz1;
441 const auto coord01_AA =
442 upper_triangle_index_ord(12, coord0_A, coord1_A);
443 const auto src = primdata_[0].targets[coord01_AA];
444 for (auto i = 0; i != ncart12; ++i) tgt[i] = src[i];
445 }
446
447 // d2 / dOi dOj += d2 / dAi dBj
448 {
449 auto coord0_A = xyz0;
450 auto coord1_B = 3 + xyz1;
451 const auto coord01_AB =
452 upper_triangle_index_ord(12, coord0_A, coord1_B);
453 const auto src = primdata_[0].targets[coord01_AB];
454 for (auto i = 0; i != ncart12; ++i) tgt[i] += src[i];
455 }
456
457 // d2 / dOi dOj += d2 / dBi dAj
458 {
459 auto coord0_B = 3 + xyz0;
460 auto coord1_A = xyz1;
461 const auto coord01_BA =
462 upper_triangle_index_ord(12, coord1_A, coord0_B);
463 const auto src = primdata_[0].targets[coord01_BA];
464 for (auto i = 0; i != ncart12; ++i) tgt[i] += src[i];
465 }
466
467 // d2 / dOi dOj += d2 / dBi dBj
468 {
469 auto coord0_B = 3 + xyz0;
470 auto coord1_B = 3 + xyz1;
471 const auto coord01_BB =
472 upper_triangle_index_ord(12, coord0_B, coord1_B);
473 const auto src = primdata_[0].targets[coord01_BB];
474 for (auto i = 0; i != ncart12; ++i) tgt[i] += src[i];
475 }
476 }
477 }
478 } // operator derivs
479
480 } break;
481
482 default: {
483 assert(deriv_order_ <= 2 && "feature not implemented");
484
485 // 1. since # of derivatives changes, remap derivatives computed
486 // by libint; targets_ will hold the "remapped" pointers to
487 // the data
488 using ShellSetDerivIterator =
489 libint2::FixedOrderedIntegerPartitionIterator<
490 std::vector<unsigned int>>;
491 ShellSetDerivIterator shellset_gaussian_diter(deriv_order_, 2);
492 ShellSetDerivIterator shellset_full_diter(deriv_order_,
493 nderivcenters_shset);
494 std::vector<unsigned int> full_deriv(3 * nderivcenters_shset, 0);
495 std::size_t s = 0;
496 while (shellset_gaussian_diter) { // loop over derivs computed
497 // by libint
498 const auto& s1s2_deriv = *shellset_gaussian_diter;
499 std::copy(std::begin(s1s2_deriv), std::end(s1s2_deriv),
500 std::begin(full_deriv));
501 const auto full_rank = ShellSetDerivIterator::rank(full_deriv);
502 targets_[full_rank] = primdata_[0].targets[s];
503 }
504 // use translational invariance to build derivatives w.r.t.
505 // operator centers
506 }
507
508 } // deriv_order_ switch
509 } // reconstruct derivatives
510 }
511 } // ltot != 0
512 } // pset (accumulation batches)
513
514 if (tform_to_solids) {
515 set_targets = false;
516 // where do spherical ints go?
517 auto* spherical_ints =
518 (accumulate_ints_in_scratch) ? scratch2_ : &scratch_[0];
519
520 // transform to solid harmonics, one shell set at a time:
521 // for each computed shell set ...
522 for (auto s = 0ul; s != num_shellsets_computed;
523 ++s, spherical_ints += n12) {
524 auto cartesian_ints = accumulate_ints_in_scratch
525 ? &scratch_[s * ncart12]
526 : primdata_[0].targets[s];
527 // transform
528 if (s1.contr[0].pure && s2.contr[0].pure) {
529 libint2::solidharmonics::tform(l1, l2, cartesian_ints, spherical_ints);
530 } else {
531 if (s1.contr[0].pure)
532 libint2::solidharmonics::tform_rows(l1, n2, cartesian_ints,
533 spherical_ints);
534 else
535 libint2::solidharmonics::tform_cols(n1, l2, cartesian_ints,
536 spherical_ints);
537 }
538 // .. and compute the destination
539 targets_[s] = spherical_ints;
540 } // loop cartesian shell set
541 } // tform to solids
542
543 if (set_targets) {
544 for (auto s = 0ul; s != num_shellsets_computed; ++s) {
545 auto cartesian_ints = accumulate_ints_in_scratch
546 ? &scratch_[s * ncart12]
547 : primdata_[0].targets[s];
548 targets_[s] = cartesian_ints;
549 }
550 }
551
552 if (cartesian_shell_normalization() == CartesianShellNormalization::uniform) {
553 std::array<std::reference_wrapper<const Shell>, 2> shells{s1, s2};
554 for (auto s = 0ul; s != num_shellsets_computed; ++s) {
555 uniform_normalize_cartesian_shells(const_cast<value_type*>(targets_[s]),
556 shells);
557 }
558 }
559
560 return targets_;
561}
562
563// generic _initializer
564__libint2_engine_inline void Engine::_initialize() {
565#define BOOST_PP_NBODYENGINE_MCR3_ncenter(product) \
566 BOOST_PP_TUPLE_ELEM(3, 1, product)
567
568#define BOOST_PP_NBODYENGINE_MCR3_default_ncenter(product) \
569 BOOST_PP_IIF(BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(3, 0, product), \
570 BOOST_PP_NBODY_OPERATOR_LAST_ONEBODY_INDEX), \
571 4, 2)
572
573#define BOOST_PP_NBODYENGINE_MCR3_NCENTER(product) \
574 BOOST_PP_IIF( \
575 BOOST_PP_NOT_EQUAL(BOOST_PP_NBODYENGINE_MCR3_ncenter(product), \
576 BOOST_PP_NBODYENGINE_MCR3_default_ncenter(product)), \
577 BOOST_PP_NBODYENGINE_MCR3_ncenter(product), BOOST_PP_EMPTY())
578
579#define BOOST_PP_NBODYENGINE_MCR3_OPER(product) \
580 BOOST_PP_LIST_AT(BOOST_PP_NBODY_OPERATOR_LIST, \
581 BOOST_PP_TUPLE_ELEM(3, 0, product))
582
583#define BOOST_PP_NBODYENGINE_MCR3_DERIV(product) \
584 BOOST_PP_IIF(BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(3, 2, product), 0), \
585 BOOST_PP_TUPLE_ELEM(3, 2, product), BOOST_PP_EMPTY())
586
587#define BOOST_PP_NBODYENGINE_MCR3_task(product) \
588 BOOST_PP_CAT(BOOST_PP_CAT(BOOST_PP_NBODYENGINE_MCR3_ncenter(product), \
589 BOOST_PP_NBODYENGINE_MCR3_OPER(product)), \
590 BOOST_PP_NBODYENGINE_MCR3_DERIV(product))
591
592#define BOOST_PP_NBODYENGINE_MCR3_TASK(product) \
593 BOOST_PP_IIF( \
594 BOOST_PP_CAT(LIBINT2_TASK_EXISTS_, \
595 BOOST_PP_NBODYENGINE_MCR3_task(product)), \
596 BOOST_PP_CAT(BOOST_PP_CAT(BOOST_PP_NBODYENGINE_MCR3_NCENTER(product), \
597 BOOST_PP_NBODYENGINE_MCR3_OPER(product)), \
598 BOOST_PP_NBODYENGINE_MCR3_DERIV(product)), \
599 default)
600
601#define BOOST_PP_NBODYENGINE_MCR3(r, product) \
602 if (static_cast<int>(oper_) == BOOST_PP_TUPLE_ELEM(3, 0, product) && \
603 static_cast<int>(rank(braket_)) == BOOST_PP_TUPLE_ELEM(3, 1, product) && \
604 deriv_order_ == BOOST_PP_TUPLE_ELEM(3, 2, product)) { \
605 hard_lmax_ = BOOST_PP_CAT(LIBINT2_MAX_AM_, \
606 BOOST_PP_NBODYENGINE_MCR3_TASK(product)) + \
607 1; \
608 hard_default_lmax_ = BOOST_PP_IF( \
609 BOOST_PP_IS_1(BOOST_PP_CAT(LIBINT2_CENTER_DEPENDENT_MAX_AM_, \
610 BOOST_PP_NBODYENGINE_MCR3_task(product))), \
611 BOOST_PP_CAT( \
612 LIBINT2_MAX_AM_, \
613 BOOST_PP_CAT(default, BOOST_PP_NBODYENGINE_MCR3_DERIV(product))) + \
614 1, \
615 std::numeric_limits<int>::max()); \
616 const auto lmax = BOOST_PP_IF( \
617 BOOST_PP_IS_1(BOOST_PP_CAT(LIBINT2_CENTER_DEPENDENT_MAX_AM_, \
618 BOOST_PP_NBODYENGINE_MCR3_task(product))), \
619 std::max(hard_lmax_, hard_default_lmax_), hard_lmax_); \
620 if (lmax_ >= lmax) { \
621 throw Engine::lmax_exceeded(BOOST_PP_STRINGIZE(BOOST_PP_NBODYENGINE_MCR3_TASK(product)), lmax, lmax_); \
622 } \
623 if (stack_size_ > 0) libint2_cleanup_default(&primdata_[0]); \
624 stack_size_ = LIBINT2_PREFIXED_NAME(BOOST_PP_CAT( \
625 libint2_need_memory_, BOOST_PP_NBODYENGINE_MCR3_TASK(product)))( \
626 lmax_); \
627 LIBINT2_PREFIXED_NAME( \
628 BOOST_PP_CAT(libint2_init_, BOOST_PP_NBODYENGINE_MCR3_TASK(product))) \
629 (&primdata_[0], lmax_, 0); \
630 BOOST_PP_IF(BOOST_PP_IS_1(LIBINT2_FLOP_COUNT), \
631 LIBINT2_PREFIXED_NAME(libint2_init_flopcounter)( \
632 &primdata_[0], primdata_.size()), \
633 BOOST_PP_EMPTY()); \
634 buildfnptrs_ = to_ptr1(LIBINT2_PREFIXED_NAME(BOOST_PP_CAT( \
635 libint2_build_, BOOST_PP_NBODYENGINE_MCR3_TASK(product)))); \
636 reset_scratch(); \
637 return; \
638 }
639
640 BOOST_PP_LIST_FOR_EACH_PRODUCT(
641 BOOST_PP_NBODYENGINE_MCR3, 3,
642 (BOOST_PP_NBODY_OPERATOR_INDEX_LIST, BOOST_PP_NBODY_BRAKET_RANK_LIST,
643 BOOST_PP_NBODY_DERIV_ORDER_LIST))
644
645 assert(false &&
646 "missing case in switch"); // either deriv_order_ or oper_ is wrong
647 abort();
648} // _initialize<R>()
649
650__libint2_engine_inline void Engine::initialize(size_t max_nprim) {
651 assert(libint2::initialized() && "libint is not initialized");
652 assert(deriv_order_ <= LIBINT2_MAX_DERIV_ORDER &&
653 "exceeded the max derivative order of the library");
654
655 // validate braket
656#ifndef INCLUDE_ONEBODY
657 assert(braket_ != BraKet::x_x &&
658 "this braket type not supported by the library; give --enable-1body "
659 "to configure");
660#endif
661#ifndef INCLUDE_ERI
662 assert(braket_ != BraKet::xx_xx &&
663 "this braket type not supported by the library; give --enable-eri to "
664 "configure");
665#endif
666#ifndef INCLUDE_ERI3
667 assert((braket_ != BraKet::xs_xx && braket_ != BraKet::xx_xs) &&
668 "this braket type not supported by the library; give --enable-eri3 to "
669 "configure");
670#endif
671#ifndef INCLUDE_ERI2
672 assert(braket_ != BraKet::xs_xs &&
673 "this braket type not supported by the library; give --enable-eri2 to "
674 "configure");
675#endif
676
677 // make sure it's no default initialized
678 if (lmax_ < 0) throw using_default_initialized();
679
680 // initialize braket, if needed
681 if (braket_ == BraKet::invalid) braket_ = default_braket(oper_);
682
683 if (max_nprim != 0) primdata_.resize(std::pow(max_nprim, braket_rank()));
684
685 // initialize targets
686 {
687 decltype(targets_)::allocator_type alloc(primdata_[0].targets);
688 targets_ = decltype(targets_)(alloc);
689 // in some cases extra memory use can be avoided if targets_ manages its own
690 // memory
691 // the only instance is where we permute derivative integrals, this calls
692 // for permuting
693 // target indices.
694 const auto permutable_targets =
695 deriv_order_ > 0 &&
696 (braket_ == BraKet::xx_xx || braket_ == BraKet::xs_xx ||
697 braket_ == BraKet::xx_xs);
698 if (permutable_targets)
699 targets_.reserve(max_ntargets + 1);
700 else
701 targets_.reserve(max_ntargets);
702 // will be resized to appropriate size in reset_scratch via _initialize
703 }
704
705#ifdef LIBINT2_ENGINE_TIMERS
706 timers.set_now_overhead(25);
707#endif
708#ifdef LIBINT2_PROFILE
709 primdata_[0].timers->set_now_overhead(25);
710#endif
711
712 _initialize();
713}
714
715namespace detail {
716__libint2_engine_inline std::vector<Engine::compute2_ptr_type>
717init_compute2_ptrs() {
718 auto max_ncompute2_ptrs = nopers_2body * nbrakets_2body * nderivorders_2body;
719 std::vector<Engine::compute2_ptr_type> result(max_ncompute2_ptrs, nullptr);
720
721#define BOOST_PP_NBODYENGINE_MCR7(r, product) \
722 if (BOOST_PP_TUPLE_ELEM(3, 0, product) >= \
723 static_cast<int>(Operator::first_2body_oper) && \
724 BOOST_PP_TUPLE_ELEM(3, 0, product) <= \
725 static_cast<int>(Operator::last_2body_oper) && \
726 BOOST_PP_TUPLE_ELEM(3, 1, product) >= \
727 static_cast<int>(BraKet::first_2body_braket) && \
728 BOOST_PP_TUPLE_ELEM(3, 1, product) <= \
729 static_cast<int>(BraKet::last_2body_braket)) { \
730 auto compute_ptr_idx = ((BOOST_PP_TUPLE_ELEM(3, 0, product) - \
731 static_cast<int>(Operator::first_2body_oper)) * \
732 nbrakets_2body + \
733 (BOOST_PP_TUPLE_ELEM(3, 1, product) - \
734 static_cast<int>(BraKet::first_2body_braket))) * \
735 nderivorders_2body + \
736 BOOST_PP_TUPLE_ELEM(3, 2, product); \
737 result.at(compute_ptr_idx) = &Engine::compute2< \
738 static_cast<Operator>(BOOST_PP_TUPLE_ELEM(3, 0, product)), \
739 static_cast<BraKet>(BOOST_PP_TUPLE_ELEM(3, 1, product)), \
740 BOOST_PP_TUPLE_ELEM(3, 2, product)>; \
741 }
742
743 BOOST_PP_LIST_FOR_EACH_PRODUCT(
744 BOOST_PP_NBODYENGINE_MCR7, 3,
745 (BOOST_PP_NBODY_OPERATOR_INDEX_LIST, BOOST_PP_NBODY_BRAKET_INDEX_LIST,
746 BOOST_PP_NBODY_DERIV_ORDER_LIST))
747
748 return result;
749}
750} // namespace detail
751
752__libint2_engine_inline const std::vector<Engine::compute2_ptr_type>&
753Engine::compute2_ptrs() const {
754 static std::vector<compute2_ptr_type> compute2_ptrs_ =
755 detail::init_compute2_ptrs();
756 return compute2_ptrs_;
757}
758
759__libint2_engine_inline unsigned int Engine::nparams() const {
760 switch (oper_) {
761 case Operator::nuclear:
762 case Operator::opVop:
763 return any_cast<
764 const operator_traits<Operator::nuclear>::oper_params_type&>(
765 params_)
766 .size();
767 case Operator::erf_nuclear:
768 case Operator::erfc_nuclear:
769 return std::get<1>(
770 any_cast<const operator_traits<
771 Operator::erfc_nuclear>::oper_params_type&>(params_))
772 .size();
773 default:
774 return 1;
775 }
776 return 1;
777}
778__libint2_engine_inline unsigned int Engine::nopers() const {
779 switch (static_cast<int>(oper_)) {
780#define BOOST_PP_NBODYENGINE_MCR4(r, data, i, elem) \
781 case i: \
782 return operator_traits<static_cast<Operator>(i)>::nopers;
783 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR4, _,
784 BOOST_PP_NBODY_OPERATOR_LIST)
785 default:
786 break;
787 }
788 assert(false && "missing case in switch"); // unreachable
789 abort();
790}
791__libint2_engine_inline unsigned int Engine::intrinsic_deriv_order() const {
792 switch (static_cast<int>(oper_)) {
793#define BOOST_PP_NBODYENGINE_MCR9(r, data, i, elem) \
794 case i: \
795 return operator_traits<static_cast<Operator>(i)>::intrinsic_deriv_order;
796 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR9, _,
797 BOOST_PP_NBODY_OPERATOR_LIST)
798 default:
799 break;
800 }
801 assert(false && "missing case in switch"); // unreachable
802 abort();
803}
804
805template <>
806__libint2_engine_inline any Engine::enforce_params_type<any>(
807 Operator oper, const any& params, bool throw_if_wrong_type) {
808 any result;
809 switch (static_cast<int>(oper)) {
810#define BOOST_PP_NBODYENGINE_MCR5A(r, data, i, elem) \
811 case i: \
812 if (any_cast<operator_traits<static_cast<Operator>(i)>::oper_params_type>( \
813 &params) != nullptr) { \
814 result = params; \
815 } else { \
816 if (throw_if_wrong_type) throw bad_any_cast(); \
817 result = operator_traits<static_cast<Operator>(i)>::default_params(); \
818 } \
819 break;
820
821 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR5A, _,
822 BOOST_PP_NBODY_OPERATOR_LIST)
823
824 default:
825 assert(false && "missing case in switch"); // missed a case?
826 abort();
827 }
828 return result;
829}
830
831template <typename Params>
832__libint2_engine_inline any Engine::enforce_params_type(
833 Operator oper, const Params& params, bool throw_if_wrong_type) {
834 any result;
835 switch (static_cast<int>(oper)) {
836#define BOOST_PP_NBODYENGINE_MCR5B(r, data, i, elem) \
837 case i: \
838 if (std::is_same<Params, operator_traits<static_cast<Operator>( \
839 i)>::oper_params_type>::value) { \
840 result = params; \
841 } else { \
842 if (throw_if_wrong_type) throw std::bad_cast(); \
843 result = operator_traits<static_cast<Operator>(i)>::default_params(); \
844 } \
845 break;
846
847 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR5B, _,
848 BOOST_PP_NBODY_OPERATOR_LIST)
849
850 default:
851 assert(false && "missing case in switch"); // missed a case?
852 abort();
853 }
854 return result;
855}
856
861__libint2_engine_inline any Engine::make_core_eval_pack(Operator oper) const {
862 any result;
863 switch (static_cast<int>(oper)) {
864#define BOOST_PP_NBODYENGINE_MCR6(r, data, i, elem) \
865 case i: \
866 result = libint2::detail::make_compressed_pair( \
867 operator_traits<static_cast<Operator>(i)>::core_eval_type::instance( \
868 braket_rank() * lmax_ + deriv_order_ + \
869 operator_traits<static_cast<Operator>( \
870 i)>::intrinsic_deriv_order, \
871 std::numeric_limits<scalar_type>::epsilon()), \
872 libint2::detail::CoreEvalScratch< \
873 operator_traits<static_cast<Operator>(i)>::core_eval_type>( \
874 braket_rank() * lmax_ + deriv_order_ + \
875 operator_traits<static_cast<Operator>( \
876 i)>::intrinsic_deriv_order)); \
877 assert(any_cast<detail::core_eval_pack_type<static_cast<Operator>(i)>>( \
878 &result) != nullptr); \
879 break;
880
881 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR6, _,
882 BOOST_PP_NBODY_OPERATOR_LIST)
883
884 default:
885 assert(false && "missing case in switch"); // missed a case?
886 abort();
887 }
888 return result;
889}
890
891__libint2_engine_inline void Engine::init_core_ints_params(const any& params) {
892 if (oper_ == Operator::delcgtg2) {
893 // [g12,[- \Del^2, g12] = 2 (\Del g12) \cdot (\Del g12)
894 // (\Del exp(-a r_12^2) \cdot (\Del exp(-b r_12^2) = 4 a b (r_{12}^2 exp(-
895 // (a+b) r_{12}^2) )
896 // i.e. need to scale each coefficient by 4 a b
897 const auto& oparams =
898 any_cast<const operator_traits<Operator::delcgtg2>::oper_params_type&>(
899 params);
900 const auto ng = oparams.size();
901 operator_traits<Operator::delcgtg2>::oper_params_type core_ints_params;
902 core_ints_params.reserve(ng * (ng + 1) / 2);
903 for (size_t b = 0; b < ng; ++b)
904 for (size_t k = 0; k <= b; ++k) {
905 const auto gexp = oparams[b].first + oparams[k].first;
906 const auto gcoeff = oparams[b].second * oparams[k].second *
907 (b == k ? 1 : 2); // if a != b include ab and ba
908 const auto gcoeff_rescaled =
909 4 * oparams[b].first * oparams[k].first * gcoeff;
910 core_ints_params.push_back(std::make_pair(gexp, gcoeff_rescaled));
911 }
912 core_ints_params_ = core_ints_params;
913 } else {
914 core_ints_params_ = params;
915 }
916}
917
918__libint2_engine_inline void Engine::compute_primdata(Libint_t& primdata,
919 const Shell& s1,
920 const Shell& s2,
921 size_t p1, size_t p2,
922 size_t oset) {
923 const auto& A = s1.O;
924 const auto& B = s2.O;
925
926 const auto alpha1 = s1.alpha[p1];
927 const auto alpha2 = s2.alpha[p2];
928
929 const auto c1 = s1.contr[0].coeff[p1];
930 const auto c2 = s2.contr[0].coeff[p2];
931
932 const auto gammap = alpha1 + alpha2;
933 const auto oogammap = 1 / gammap;
934 const auto rhop_over_alpha1 = alpha2 * oogammap;
935 const auto rhop = alpha1 * rhop_over_alpha1;
936 const auto Px = (alpha1 * A[0] + alpha2 * B[0]) * oogammap;
937 const auto Py = (alpha1 * A[1] + alpha2 * B[1]) * oogammap;
938 const auto Pz = (alpha1 * A[2] + alpha2 * B[2]) * oogammap;
939 const auto AB_x = A[0] - B[0];
940 const auto AB_y = A[1] - B[1];
941 const auto AB_z = A[2] - B[2];
942 const auto AB2_x = AB_x * AB_x;
943 const auto AB2_y = AB_y * AB_y;
944 const auto AB2_z = AB_z * AB_z;
945
946 assert(LIBINT2_SHELLQUARTET_SET == LIBINT2_SHELLQUARTET_SET_STANDARD &&
947 "non-standard shell ordering");
948
949 const auto oper_is_nuclear =
950 (oper_ == Operator::nuclear || oper_ == Operator::erf_nuclear ||
951 oper_ == Operator::erfc_nuclear || oper_ == Operator::opVop);
952
953 // need to use HRR? see strategy.cc
954 const auto l1 = s1.contr[0].l;
955 const auto l2 = s2.contr[0].l;
956 const bool use_hrr =
957 (oper_is_nuclear || oper_ == Operator::sphemultipole) && l1 > 0 && l2 > 0;
958 // unlike the 2-body ints, can go both ways, determine which way to go (the
959 // logic must match TwoCenter_OS_Tactic)
960 const bool hrr_ket_to_bra = l1 >= l2;
961 if (use_hrr) {
962 if (hrr_ket_to_bra) {
963#if LIBINT2_DEFINED(eri, AB_x)
964 primdata.AB_x[0] = AB_x;
965#endif
966#if LIBINT2_DEFINED(eri, AB_y)
967 primdata.AB_y[0] = AB_y;
968#endif
969#if LIBINT2_DEFINED(eri, AB_z)
970 primdata.AB_z[0] = AB_z;
971#endif
972 } else {
973#if LIBINT2_DEFINED(eri, BA_x)
974 primdata.BA_x[0] = -AB_x;
975#endif
976#if LIBINT2_DEFINED(eri, BA_y)
977 primdata.BA_y[0] = -AB_y;
978#endif
979#if LIBINT2_DEFINED(eri, BA_z)
980 primdata.BA_z[0] = -AB_z;
981#endif
982 }
983 }
984
985 // figure out whether will do VRR on center A and/or B
986// if ((!use_hrr && l1 > 0) || hrr_ket_to_bra) {
987#if LIBINT2_DEFINED(eri, PA_x)
988 primdata.PA_x[0] = Px - A[0];
989#endif
990#if LIBINT2_DEFINED(eri, PA_y)
991 primdata.PA_y[0] = Py - A[1];
992#endif
993#if LIBINT2_DEFINED(eri, PA_z)
994 primdata.PA_z[0] = Pz - A[2];
995#endif
996// }
997//
998// if ((!use_hrr && l2 > 0) || !hrr_ket_to_bra) {
999#if LIBINT2_DEFINED(eri, PB_x)
1000 primdata.PB_x[0] = Px - B[0];
1001#endif
1002#if LIBINT2_DEFINED(eri, PB_y)
1003 primdata.PB_y[0] = Py - B[1];
1004#endif
1005#if LIBINT2_DEFINED(eri, PB_z)
1006 primdata.PB_z[0] = Pz - B[2];
1007#endif
1008 // }
1009
1010 if (oper_ == Operator::emultipole1 || oper_ == Operator::emultipole2 ||
1011 oper_ == Operator::emultipole3) {
1012 const auto& O = any_cast<
1013 const operator_traits<Operator::emultipole1>::oper_params_type&>(
1014 params_); // same as emultipoleX
1015#if LIBINT2_DEFINED(eri, BO_x)
1016 primdata.BO_x[0] = B[0] - O[0];
1017#endif
1018#if LIBINT2_DEFINED(eri, BO_y)
1019 primdata.BO_y[0] = B[1] - O[1];
1020#endif
1021#if LIBINT2_DEFINED(eri, BO_z)
1022 primdata.BO_z[0] = B[2] - O[2];
1023#endif
1024 }
1025 if (oper_ == Operator::sphemultipole) {
1026 const auto& O = any_cast<
1027 const operator_traits<Operator::emultipole1>::oper_params_type&>(
1028 params_);
1029#if LIBINT2_DEFINED(eri, PO_x)
1030 primdata.PO_x[0] = Px - O[0];
1031#endif
1032#if LIBINT2_DEFINED(eri, PO_y)
1033 primdata.PO_y[0] = Py - O[1];
1034#endif
1035#if LIBINT2_DEFINED(eri, PO_z)
1036 primdata.PO_z[0] = Pz - O[2];
1037#endif
1038#if LIBINT2_DEFINED(eri, PO2)
1039 primdata.PO2[0] = (Px - O[0]) * (Px - O[0]) + (Py - O[1]) * (Py - O[1]) +
1040 (Pz - O[2]) * (Pz - O[2]);
1041#endif
1042 }
1043
1044#if LIBINT2_DEFINED(eri, oo2z)
1045 primdata.oo2z[0] = 0.5 * oogammap;
1046#endif
1047
1048 decltype(c1) sqrt_PI(1.77245385090551602729816748334);
1049 const auto xyz_pfac = sqrt_PI * sqrt(oogammap);
1050 const auto ovlp_ss_x = exp(-rhop * AB2_x) * xyz_pfac * c1 * c2 * scale_;
1051 const auto ovlp_ss_y = exp(-rhop * AB2_y) * xyz_pfac;
1052 const auto ovlp_ss_z = exp(-rhop * AB2_z) * xyz_pfac;
1053
1054 primdata._0_Overlap_0_x[0] = ovlp_ss_x;
1055 primdata._0_Overlap_0_y[0] = ovlp_ss_y;
1056 primdata._0_Overlap_0_z[0] = ovlp_ss_z;
1057
1058 if (oper_ == Operator::kinetic || (deriv_order_ > 0) ||
1059 oper_ == Operator::opVop) {
1060#if LIBINT2_DEFINED(eri, two_alpha0_bra)
1061 primdata.two_alpha0_bra[0] = 2.0 * alpha1;
1062#endif
1063#if LIBINT2_DEFINED(eri, two_alpha0_ket)
1064 primdata.two_alpha0_ket[0] = 2.0 * alpha2;
1065#endif
1066 }
1067
1068 if (oper_is_nuclear) {
1069 const auto& params =
1070 (oper_ == Operator::nuclear || oper_ == Operator::opVop)
1071 ? any_cast<
1072 const operator_traits<Operator::nuclear>::oper_params_type&>(
1073 params_)
1074 : std::get<1>(
1075 any_cast<const operator_traits<
1076 Operator::erfc_nuclear>::oper_params_type&>(params_));
1077
1078 const auto& C = params[oset].second;
1079 const auto& q = params[oset].first;
1080#if LIBINT2_DEFINED(eri, PC_x)
1081 primdata.PC_x[0] = Px - C[0];
1082#endif
1083#if LIBINT2_DEFINED(eri, PC_y)
1084 primdata.PC_y[0] = Py - C[1];
1085#endif
1086#if LIBINT2_DEFINED(eri, PC_z)
1087 primdata.PC_z[0] = Pz - C[2];
1088#endif
1089
1090#if LIBINT2_DEFINED(eri, rho12_over_alpha1) || \
1091 LIBINT2_DEFINED(eri, rho12_over_alpha2)
1092 if (deriv_order_ + intrinsic_deriv_order() > 0) {
1093#if LIBINT2_DEFINED(eri, rho12_over_alpha1)
1094 primdata.rho12_over_alpha1[0] = rhop_over_alpha1;
1095#endif
1096#if LIBINT2_DEFINED(eri, rho12_over_alpha2)
1097 primdata.rho12_over_alpha2[0] = alpha1 * oogammap;
1098#endif
1099 }
1100#endif
1101#if LIBINT2_DEFINED(eri, PC_x) && LIBINT2_DEFINED(eri, PC_y) && \
1102 LIBINT2_DEFINED(eri, PC_z)
1103 const auto PC2 = primdata.PC_x[0] * primdata.PC_x[0] +
1104 primdata.PC_y[0] * primdata.PC_y[0] +
1105 primdata.PC_z[0] * primdata.PC_z[0];
1106 const scalar_type U = gammap * PC2;
1107 const auto mmax =
1108 s1.contr[0].l + s2.contr[0].l + deriv_order_ + intrinsic_deriv_order();
1109 auto* fm_ptr = &(primdata.LIBINT_T_S_ELECPOT_S(0)[0]);
1110 if (oper_ == Operator::nuclear || oper_ == Operator::opVop) {
1111 const auto& fm_engine_ptr =
1112 any_cast<const detail::core_eval_pack_type<Operator::nuclear>&>(
1113 core_eval_pack_)
1114 .first();
1115 fm_engine_ptr->eval(fm_ptr, U, mmax);
1116 } else if (oper_ == Operator::erf_nuclear) {
1117 const auto& core_eval_ptr =
1118 any_cast<const detail::core_eval_pack_type<Operator::erf_nuclear>&>(
1119 core_eval_pack_)
1120 .first();
1121 const auto& core_ints_params = std::get<0>(
1122 any_cast<const typename operator_traits<
1123 Operator::erf_nuclear>::oper_params_type&>(core_ints_params_));
1124 core_eval_ptr->eval(fm_ptr, gammap, U, mmax, core_ints_params);
1125 } else if (oper_ == Operator::erfc_nuclear) {
1126 const auto& core_eval_ptr =
1127 any_cast<const detail::core_eval_pack_type<Operator::erfc_nuclear>&>(
1128 core_eval_pack_)
1129 .first();
1130 const auto& core_ints_params = std::get<0>(
1131 any_cast<const typename operator_traits<
1132 Operator::erfc_nuclear>::oper_params_type&>(core_ints_params_));
1133 core_eval_ptr->eval(fm_ptr, gammap, U, mmax, core_ints_params);
1134 }
1135
1136 decltype(U) two_o_sqrt_PI(1.12837916709551257389615890312);
1137 const auto pfac =
1138 -q * sqrt(gammap) * two_o_sqrt_PI * ovlp_ss_x * ovlp_ss_y * ovlp_ss_z;
1139 const auto m_fence = mmax + 1;
1140 for (auto m = 0; m != m_fence; ++m) {
1141 fm_ptr[m] *= pfac;
1142 }
1143#endif
1144 }
1145} // Engine::compute_primdata()
1146
1151template <Operator op, BraKet bk, size_t der>
1152__libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2(
1153 const libint2::Shell& tbra1, const libint2::Shell& tbra2,
1154 const libint2::Shell& tket1, const libint2::Shell& tket2,
1155 const ShellPair* tspbra, const ShellPair* tspket) {
1156 assert(op == oper_ && "Engine::compute2 -- operator mismatch");
1157 assert(bk == braket_ && "Engine::compute2 -- braket mismatch");
1158 assert(der == deriv_order_ && "Engine::compute2 -- deriv_order mismatch");
1159 assert(((tspbra == nullptr && tspket == nullptr) ||
1160 (tspbra != nullptr && tspket != nullptr)) &&
1161 "Engine::compute2 -- expects zero or two ShellPair objects");
1162 assert(screening_method_ != ScreeningMethod::Invalid);
1163
1164 //
1165 // i.e. bra and ket refer to chemists bra and ket
1166 //
1167
1168 // can only handle 1 contraction at a time
1169 assert((tbra1.ncontr() == 1 && tbra2.ncontr() == 1 && tket1.ncontr() == 1 &&
1170 tket2.ncontr() == 1) &&
1171 "generally-contracted shells are not yet supported");
1172
1173 // angular momentum limit obeyed? can only be fully checked in
1174 // braket-dependent code, here only do a basic test that does not guarantee
1175 // that the shell-set can be computed
1176 assert(tbra1.contr[0].l <= lmax_ && "the angular momentum limit is exceeded");
1177 assert(tbra2.contr[0].l <= lmax_ && "the angular momentum limit is exceeded");
1178 assert(tket1.contr[0].l <= lmax_ && "the angular momentum limit is exceeded");
1179 assert(tket2.contr[0].l <= lmax_ && "the angular momentum limit is exceeded");
1180
1181#if LIBINT2_SHELLQUARTET_SET == \
1182 LIBINT2_SHELLQUARTET_SET_STANDARD // standard angular momentum ordering
1183 const auto swap_tbra = (tbra1.contr[0].l < tbra2.contr[0].l);
1184 const auto swap_tket = (tket1.contr[0].l < tket2.contr[0].l);
1185 const auto swap_braket =
1186 ((braket_ == BraKet::xx_xx) && (tbra1.contr[0].l + tbra2.contr[0].l >
1187 tket1.contr[0].l + tket2.contr[0].l)) ||
1188 braket_ == BraKet::xx_xs;
1189#else // orca angular momentum ordering
1190 const auto swap_tbra = (tbra1.contr[0].l > tbra2.contr[0].l);
1191 const auto swap_tket = (tket1.contr[0].l > tket2.contr[0].l);
1192 const auto swap_braket =
1193 ((braket_ == BraKet::xx_xx) && (tbra1.contr[0].l + tbra2.contr[0].l <
1194 tket1.contr[0].l + tket2.contr[0].l)) ||
1195 braket_ == BraKet::xx_xs;
1196 assert(false && "feature not implemented");
1197 abort();
1198#endif
1199 const auto& bra1 =
1200 swap_braket ? (swap_tket ? tket2 : tket1) : (swap_tbra ? tbra2 : tbra1);
1201 const auto& bra2 =
1202 swap_braket ? (swap_tket ? tket1 : tket2) : (swap_tbra ? tbra1 : tbra2);
1203 const auto& ket1 =
1204 swap_braket ? (swap_tbra ? tbra2 : tbra1) : (swap_tket ? tket2 : tket1);
1205 const auto& ket2 =
1206 swap_braket ? (swap_tbra ? tbra1 : tbra2) : (swap_tket ? tket1 : tket2);
1207 const auto swap_bra = swap_braket ? swap_tket : swap_tbra;
1208 const auto swap_ket = swap_braket ? swap_tbra : swap_tket;
1209 // "permute" also the user-provided shell pair data
1210 const auto* spbra_precomputed = swap_braket ? tspket : tspbra;
1211 const auto* spket_precomputed = swap_braket ? tspbra : tspket;
1212 assert(((spbra_precomputed && spket_precomputed) ||
1213 (screening_method_ == ScreeningMethod::Original ||
1214 screening_method_ == ScreeningMethod::Conservative)) &&
1215 "Engine::compute2: without precomputed shell pair data can only use "
1216 "original or conservative screening methods");
1217
1218 const auto tform = bra1.contr[0].pure || bra2.contr[0].pure ||
1219 ket1.contr[0].pure || ket2.contr[0].pure;
1220 const auto permute = swap_braket || swap_tbra || swap_tket;
1221 const auto use_scratch = permute || tform;
1222
1223 // assert # of primitive pairs
1224 auto nprim_bra1 = bra1.nprim();
1225 auto nprim_bra2 = bra2.nprim();
1226 auto nprim_ket1 = ket1.nprim();
1227 auto nprim_ket2 = ket2.nprim();
1228
1229 // adjust max angular momentum, if needed
1230 auto lmax = std::max(std::max(bra1.contr[0].l, bra2.contr[0].l),
1231 std::max(ket1.contr[0].l, ket2.contr[0].l));
1232 assert(lmax <= lmax_ && "the angular momentum limit is exceeded");
1233 const auto lmax_bra = std::max(bra1.contr[0].l, bra2.contr[0].l);
1234 const auto lmax_ket = std::max(ket1.contr[0].l, ket2.contr[0].l);
1235
1236#ifdef LIBINT2_ENGINE_PROFILE_CLASS
1237 class_id id(bra1.contr[0].l, bra2.contr[0].l, ket1.contr[0].l,
1238 ket2.contr[0].l);
1239 if (class_profiles.find(id) == class_profiles.end()) {
1240 class_profile dummy;
1241 class_profiles[id] = dummy;
1242 }
1243#endif
1244
1245// compute primitive data
1246#ifdef LIBINT2_ENGINE_TIMERS
1247 timers.start(0);
1248#endif
1249 {
1250 auto p = 0;
1251 // initialize shell pairs, if not given ...
1252 // since screening primitive pairs for bra is not possible without knowing
1253 // the worst-case primitive data in ket (and vice versa) using ln_precision_
1254 // is not safe, but should work fine for moderate basis sets ... precompute
1255 // shell pair data yourself to guarantee proper screening of primitive pairs
1256 // (see Engine::set_precision() for details of how to screen primitives for
1257 // a given target precision of the integrals)
1258 const auto target_shellpair_ln_precision = ln_precision_;
1259 const auto recompute_spbra =
1260 !spbra_precomputed ||
1261 spbra_precomputed->ln_prec > target_shellpair_ln_precision;
1262 const auto recompute_spket =
1263 !spket_precomputed ||
1264 spket_precomputed->ln_prec > target_shellpair_ln_precision;
1265 const ShellPair& spbra =
1266 recompute_spbra
1267 ? (spbra_.init(bra1, bra2, target_shellpair_ln_precision,
1268 screening_method_),
1269 spbra_)
1270 : *spbra_precomputed;
1271 const ShellPair& spket =
1272 recompute_spket
1273 ? (spket_.init(ket1, ket2, target_shellpair_ln_precision,
1274 screening_method_),
1275 spket_)
1276 : *spket_precomputed;
1277 assert(spbra.screening_method_ == screening_method_ &&
1278 spket.screening_method_ == screening_method_ &&
1279 "Engine::compute2: received ShellPair initialized for an "
1280 "incompatible screening method");
1281 // determine whether shell pair data refers to the actual ({bra1,bra2}) or
1282 // swapped ({bra2,bra1}) pairs if computed the shell pair data here then
1283 // it's always in actual order, otherwise check swap_bra/swap_ket
1284 const auto spbra_is_swapped = recompute_spbra ? false : swap_bra;
1285 const auto spket_is_swapped = recompute_spket ? false : swap_ket;
1286
1287 using real_t = Shell::real_t;
1288 // swapping bra turns AB into BA = -AB
1289 real_t BA[3];
1290 if (spbra_is_swapped) {
1291 for (auto xyz = 0; xyz != 3; ++xyz) BA[xyz] = -spbra_precomputed->AB[xyz];
1292 }
1293 const auto& AB = spbra_is_swapped ? BA : spbra.AB;
1294 // swapping ket turns CD into DC = -CD
1295 real_t DC[3];
1296 if (spket_is_swapped) {
1297 for (auto xyz = 0; xyz != 3; ++xyz) DC[xyz] = -spket_precomputed->AB[xyz];
1298 }
1299 const auto& CD = spket_is_swapped ? DC : spket.AB;
1300
1301 const auto& A = bra1.O;
1302 const auto& B = bra2.O;
1303 const auto& C = ket1.O;
1304 const auto& D = ket2.O;
1305
1306 // compute all primitive quartet data
1307 const auto npbra = spbra.primpairs.size();
1308 const auto npket = spket.primpairs.size();
1309 const scalar_type npbraket = npbra * npket;
1310 for (auto pb = 0; pb != npbra; ++pb) {
1311 for (auto pk = 0; pk != npket; ++pk) {
1312 // primitive quartet coarse screening:
1313 if (spbra.primpairs[pb].ln_scr + spket.primpairs[pk].ln_scr >
1314 ln_precision_) {
1315 Libint_t& primdata = primdata_[p];
1316 const auto& sbra1 = bra1;
1317 const auto& sbra2 = bra2;
1318 const auto& sket1 = ket1;
1319 const auto& sket2 = ket2;
1320 auto pbra = pb;
1321 auto pket = pk;
1322
1323 const auto& spbrapp = spbra.primpairs[pbra];
1324 const auto& spketpp = spket.primpairs[pket];
1325 // if shell-pair data given by user
1326 const auto& pbra1 = spbra_is_swapped ? spbrapp.p2 : spbrapp.p1;
1327 const auto& pbra2 = spbra_is_swapped ? spbrapp.p1 : spbrapp.p2;
1328 const auto& pket1 = spket_is_swapped ? spketpp.p2 : spketpp.p1;
1329 const auto& pket2 = spket_is_swapped ? spketpp.p1 : spketpp.p2;
1330
1331 const auto alpha0 = sbra1.alpha[pbra1];
1332 const auto alpha1 = sbra2.alpha[pbra2];
1333 const auto alpha2 = sket1.alpha[pket1];
1334 const auto alpha3 = sket2.alpha[pket2];
1335
1336 const auto c0 = sbra1.contr[0].coeff[pbra1];
1337 const auto c1 = sbra2.contr[0].coeff[pbra2];
1338 const auto c2 = sket1.contr[0].coeff[pket1];
1339 const auto c3 = sket2.contr[0].coeff[pket2];
1340
1341 const auto l0 = sbra1.contr[0].l;
1342 const auto l1 = sbra2.contr[0].l;
1343 const auto l2 = sket1.contr[0].l;
1344 const auto l3 = sket2.contr[0].l;
1345 const auto l = l0 + l1 + l2 + l3;
1346
1347 const auto gammap = alpha0 + alpha1;
1348 const auto oogammap = spbrapp.one_over_gamma;
1349 const auto rhop = alpha0 * alpha1 * oogammap;
1350
1351 const auto gammaq = alpha2 + alpha3;
1352 const auto oogammaq = spketpp.one_over_gamma;
1353 const auto rhoq = alpha2 * alpha3 * oogammaq;
1354
1355 const auto& P = spbrapp.P;
1356 const auto& Q = spketpp.P;
1357 const auto PQx = P[0] - Q[0];
1358 const auto PQy = P[1] - Q[1];
1359 const auto PQz = P[2] - Q[2];
1360 const auto PQ2 = PQx * PQx + PQy * PQy + PQz * PQz;
1361
1362 const auto K12 = spbrapp.K * spketpp.K;
1363 const auto gammapq = gammap + gammaq;
1364 const auto sqrt_gammapq = sqrt(gammapq);
1365 const auto oogammapq = 1.0 / (gammapq);
1366 auto pfac = K12 * sqrt_gammapq * oogammapq;
1367 pfac *= c0 * c1 * c2 * c3 * scale_;
1368
1369 // original and conservative methods: screen primitive integral using
1370 // actual pfac
1371 if (static_cast<int>(screening_method_) &
1372 (static_cast<int>(ScreeningMethod::Original) |
1373 static_cast<int>(ScreeningMethod::Conservative))) {
1374 scalar_type magnitude_estimate = std::abs(pfac);
1375 if (screening_method_ == ScreeningMethod::Conservative) {
1376 // magnitude of primitive (ab|cd) integral for nonzero L differs
1377 // from that of (00|00) geometric and exponent-dependent factor,
1378 // some of which only depend on bra or ket, and some are bra-ket
1379 // dependent ... here we account only for bra-only and ket-only
1380 // nonspherical factors
1381 const auto nonspherical_pfac_magnitude = std::max(
1382 1., spbrapp.nonsph_screen_fac * spketpp.nonsph_screen_fac);
1383 magnitude_estimate *= nonspherical_pfac_magnitude * npbraket;
1384 }
1385 if (magnitude_estimate < precision_) continue;
1386 }
1387
1388 {
1389 const scalar_type rho = gammap * gammaq * oogammapq;
1390 const scalar_type T = PQ2 * rho;
1391 auto* gm_ptr = &(primdata.LIBINT_T_SS_EREP_SS(0)[0]);
1392 const auto mmax = l + deriv_order_;
1393
1394 if (!skip_core_ints) {
1395 switch (oper_) {
1396 case Operator::coulomb: {
1397 const auto& core_eval_ptr =
1398 any_cast<const detail::core_eval_pack_type<
1399 Operator::coulomb>&>(core_eval_pack_)
1400 .first();
1401 core_eval_ptr->eval(gm_ptr, T, mmax);
1402 } break;
1403 case Operator::cgtg_x_coulomb: {
1404 const auto& core_eval_ptr =
1405 any_cast<const detail::core_eval_pack_type<
1406 Operator::cgtg_x_coulomb>&>(core_eval_pack_)
1407 .first();
1408 auto& core_eval_scratch =
1409 any_cast<detail::core_eval_pack_type<
1410 Operator::cgtg_x_coulomb>&>(core_eval_pack_)
1411 .second();
1412 const auto& core_ints_params =
1413 any_cast<const typename operator_traits<
1414 Operator::cgtg>::oper_params_type&>(
1415 core_ints_params_);
1416 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params,
1417 &core_eval_scratch);
1418 } break;
1419 case Operator::cgtg: {
1420 const auto& core_eval_ptr =
1421 any_cast<
1422 const detail::core_eval_pack_type<Operator::cgtg>&>(
1423 core_eval_pack_)
1424 .first();
1425 const auto& core_ints_params =
1426 any_cast<const typename operator_traits<
1427 Operator::cgtg>::oper_params_type&>(
1428 core_ints_params_);
1429 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params);
1430 } break;
1431 case Operator::delcgtg2: {
1432 const auto& core_eval_ptr =
1433 any_cast<const detail::core_eval_pack_type<
1434 Operator::delcgtg2>&>(core_eval_pack_)
1435 .first();
1436 const auto& core_ints_params =
1437 any_cast<const typename operator_traits<
1438 Operator::cgtg>::oper_params_type&>(
1439 core_ints_params_);
1440 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params);
1441 } break;
1442 case Operator::delta: {
1443 const auto& core_eval_ptr =
1444 any_cast<
1445 const detail::core_eval_pack_type<Operator::delta>&>(
1446 core_eval_pack_)
1447 .first();
1448 core_eval_ptr->eval(gm_ptr, rho, T, mmax);
1449 } break;
1450 case Operator::r12: {
1451 const auto& core_eval_ptr =
1452 any_cast<
1453 const detail::core_eval_pack_type<Operator::r12>&>(
1454 core_eval_pack_)
1455 .first();
1456 core_eval_ptr->eval(gm_ptr, rho, T, mmax);
1457 } break;
1458 case Operator::erf_coulomb: {
1459 const auto& core_eval_ptr =
1460 any_cast<const detail::core_eval_pack_type<
1461 Operator::erf_coulomb>&>(core_eval_pack_)
1462 .first();
1463 const auto& core_ints_params =
1464 any_cast<const typename operator_traits<
1465 Operator::erf_coulomb>::oper_params_type&>(
1466 core_ints_params_);
1467 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params);
1468 } break;
1469 case Operator::erfc_coulomb: {
1470 const auto& core_eval_ptr =
1471 any_cast<const detail::core_eval_pack_type<
1472 Operator::erfc_coulomb>&>(core_eval_pack_)
1473 .first();
1474 const auto& core_ints_params =
1475 any_cast<const typename operator_traits<
1476 Operator::erfc_coulomb>::oper_params_type&>(
1477 core_ints_params_);
1478 core_eval_ptr->eval(gm_ptr, rho, T, mmax, core_ints_params);
1479 } break;
1480 case Operator::stg: {
1481 const auto& core_eval_ptr =
1482 any_cast<
1483 const detail::core_eval_pack_type<Operator::stg>&>(
1484 core_eval_pack_)
1485 .first();
1486 auto zeta = any_cast<const typename operator_traits<
1487 Operator::stg>::oper_params_type&>(core_ints_params_);
1488 const auto one_over_rho = gammapq * oogammap * oogammaq;
1489 core_eval_ptr->eval_slater(gm_ptr, one_over_rho, T, mmax,
1490 zeta);
1491 } break;
1492 case Operator::stg_x_coulomb: {
1493 const auto& core_eval_ptr =
1494 any_cast<
1495 const detail::core_eval_pack_type<Operator::stg>&>(
1496 core_eval_pack_)
1497 .first();
1498 auto zeta = any_cast<const typename operator_traits<
1499 Operator::stg>::oper_params_type&>(core_ints_params_);
1500 const auto one_over_rho = gammapq * oogammap * oogammaq;
1501 core_eval_ptr->eval_yukawa(gm_ptr, one_over_rho, T, mmax,
1502 zeta);
1503 } break;
1504 default:
1505 assert(false && "missing case in a switch"); // unreachable
1506 abort();
1507 }
1508 }
1509
1510 for (auto m = 0; m != mmax + 1; ++m) {
1511 gm_ptr[m] *= pfac;
1512 }
1513
1514 if (mmax != 0) {
1515 if (braket_ == BraKet::xx_xx) {
1516#if LIBINT2_DEFINED(eri, PA_x)
1517 primdata.PA_x[0] = P[0] - A[0];
1518#endif
1519#if LIBINT2_DEFINED(eri, PA_y)
1520 primdata.PA_y[0] = P[1] - A[1];
1521#endif
1522#if LIBINT2_DEFINED(eri, PA_z)
1523 primdata.PA_z[0] = P[2] - A[2];
1524#endif
1525#if LIBINT2_DEFINED(eri, PB_x)
1526 primdata.PB_x[0] = P[0] - B[0];
1527#endif
1528#if LIBINT2_DEFINED(eri, PB_y)
1529 primdata.PB_y[0] = P[1] - B[1];
1530#endif
1531#if LIBINT2_DEFINED(eri, PB_z)
1532 primdata.PB_z[0] = P[2] - B[2];
1533#endif
1534 }
1535
1536 if (braket_ != BraKet::xs_xs) {
1537#if LIBINT2_DEFINED(eri, QC_x)
1538 primdata.QC_x[0] = Q[0] - C[0];
1539#endif
1540#if LIBINT2_DEFINED(eri, QC_y)
1541 primdata.QC_y[0] = Q[1] - C[1];
1542#endif
1543#if LIBINT2_DEFINED(eri, QC_z)
1544 primdata.QC_z[0] = Q[2] - C[2];
1545#endif
1546#if LIBINT2_DEFINED(eri, QD_x)
1547 primdata.QD_x[0] = Q[0] - D[0];
1548#endif
1549#if LIBINT2_DEFINED(eri, QD_y)
1550 primdata.QD_y[0] = Q[1] - D[1];
1551#endif
1552#if LIBINT2_DEFINED(eri, QD_z)
1553 primdata.QD_z[0] = Q[2] - D[2];
1554#endif
1555 }
1556
1557 if (braket_ == BraKet::xx_xx) {
1558#if LIBINT2_DEFINED(eri, AB_x)
1559 primdata.AB_x[0] = AB[0];
1560#endif
1561#if LIBINT2_DEFINED(eri, AB_y)
1562 primdata.AB_y[0] = AB[1];
1563#endif
1564#if LIBINT2_DEFINED(eri, AB_z)
1565 primdata.AB_z[0] = AB[2];
1566#endif
1567#if LIBINT2_DEFINED(eri, BA_x)
1568 primdata.BA_x[0] = -AB[0];
1569#endif
1570#if LIBINT2_DEFINED(eri, BA_y)
1571 primdata.BA_y[0] = -AB[1];
1572#endif
1573#if LIBINT2_DEFINED(eri, BA_z)
1574 primdata.BA_z[0] = -AB[2];
1575#endif
1576 }
1577
1578 if (braket_ != BraKet::xs_xs) {
1579#if LIBINT2_DEFINED(eri, CD_x)
1580 primdata.CD_x[0] = CD[0];
1581#endif
1582#if LIBINT2_DEFINED(eri, CD_y)
1583 primdata.CD_y[0] = CD[1];
1584#endif
1585#if LIBINT2_DEFINED(eri, CD_z)
1586 primdata.CD_z[0] = CD[2];
1587#endif
1588#if LIBINT2_DEFINED(eri, DC_x)
1589 primdata.DC_x[0] = -CD[0];
1590#endif
1591#if LIBINT2_DEFINED(eri, DC_y)
1592 primdata.DC_y[0] = -CD[1];
1593#endif
1594#if LIBINT2_DEFINED(eri, DC_z)
1595 primdata.DC_z[0] = -CD[2];
1596#endif
1597 }
1598
1599 const auto gammap_o_gammapgammaq = oogammapq * gammap;
1600 const auto gammaq_o_gammapgammaq = oogammapq * gammaq;
1601
1602 const auto Wx =
1603 (gammap_o_gammapgammaq * P[0] + gammaq_o_gammapgammaq * Q[0]);
1604 const auto Wy =
1605 (gammap_o_gammapgammaq * P[1] + gammaq_o_gammapgammaq * Q[1]);
1606 const auto Wz =
1607 (gammap_o_gammapgammaq * P[2] + gammaq_o_gammapgammaq * Q[2]);
1608
1609 if (deriv_order_ > 0 || lmax_bra > 0) {
1610#if LIBINT2_DEFINED(eri, WP_x)
1611 primdata.WP_x[0] = Wx - P[0];
1612#endif
1613#if LIBINT2_DEFINED(eri, WP_y)
1614 primdata.WP_y[0] = Wy - P[1];
1615#endif
1616#if LIBINT2_DEFINED(eri, WP_z)
1617 primdata.WP_z[0] = Wz - P[2];
1618#endif
1619 }
1620 if (deriv_order_ > 0 || lmax_ket > 0) {
1621#if LIBINT2_DEFINED(eri, WQ_x)
1622 primdata.WQ_x[0] = Wx - Q[0];
1623#endif
1624#if LIBINT2_DEFINED(eri, WQ_y)
1625 primdata.WQ_y[0] = Wy - Q[1];
1626#endif
1627#if LIBINT2_DEFINED(eri, WQ_z)
1628 primdata.WQ_z[0] = Wz - Q[2];
1629#endif
1630 }
1631#if LIBINT2_DEFINED(eri, oo2z)
1632 primdata.oo2z[0] = 0.5 * oogammap;
1633#endif
1634#if LIBINT2_DEFINED(eri, oo2e)
1635 primdata.oo2e[0] = 0.5 * oogammaq;
1636#endif
1637#if LIBINT2_DEFINED(eri, oo2ze)
1638 primdata.oo2ze[0] = 0.5 * oogammapq;
1639#endif
1640#if LIBINT2_DEFINED(eri, roz)
1641 primdata.roz[0] = rho * oogammap;
1642#endif
1643#if LIBINT2_DEFINED(eri, roe)
1644 primdata.roe[0] = rho * oogammaq;
1645#endif
1646
1647// using ITR?
1648#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_0_x)
1649 primdata.TwoPRepITR_pfac0_0_0_x[0] =
1650 -(alpha1 * AB[0] + alpha3 * CD[0]) * oogammap;
1651#endif
1652#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_0_y)
1653 primdata.TwoPRepITR_pfac0_0_0_y[0] =
1654 -(alpha1 * AB[1] + alpha3 * CD[1]) * oogammap;
1655#endif
1656#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_0_z)
1657 primdata.TwoPRepITR_pfac0_0_0_z[0] =
1658 -(alpha1 * AB[2] + alpha3 * CD[2]) * oogammap;
1659#endif
1660#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_0_x)
1661 primdata.TwoPRepITR_pfac0_1_0_x[0] =
1662 -(alpha1 * AB[0] + alpha3 * CD[0]) * oogammaq;
1663#endif
1664#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_0_y)
1665 primdata.TwoPRepITR_pfac0_1_0_y[0] =
1666 -(alpha1 * AB[1] + alpha3 * CD[1]) * oogammaq;
1667#endif
1668#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_0_z)
1669 primdata.TwoPRepITR_pfac0_1_0_z[0] =
1670 -(alpha1 * AB[2] + alpha3 * CD[2]) * oogammaq;
1671#endif
1672#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_1_x)
1673 primdata.TwoPRepITR_pfac0_0_1_x[0] =
1674 (alpha0 * AB[0] + alpha2 * CD[0]) * oogammap;
1675#endif
1676#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_1_y)
1677 primdata.TwoPRepITR_pfac0_0_1_y[0] =
1678 (alpha0 * AB[1] + alpha2 * CD[1]) * oogammap;
1679#endif
1680#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_0_1_z)
1681 primdata.TwoPRepITR_pfac0_0_1_z[0] =
1682 (alpha0 * AB[2] + alpha2 * CD[2]) * oogammap;
1683#endif
1684#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_1_x)
1685 primdata.TwoPRepITR_pfac0_1_1_x[0] =
1686 (alpha0 * AB[0] + alpha2 * CD[0]) * oogammaq;
1687#endif
1688#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_1_y)
1689 primdata.TwoPRepITR_pfac0_1_1_y[0] =
1690 (alpha0 * AB[1] + alpha2 * CD[1]) * oogammaq;
1691#endif
1692#if LIBINT2_DEFINED(eri, TwoPRepITR_pfac0_1_1_z)
1693 primdata.TwoPRepITR_pfac0_1_1_z[0] =
1694 (alpha0 * AB[2] + alpha2 * CD[2]) * oogammaq;
1695#endif
1696#if LIBINT2_DEFINED(eri, eoz)
1697 primdata.eoz[0] = gammaq * oogammap;
1698#endif
1699#if LIBINT2_DEFINED(eri, zoe)
1700 primdata.zoe[0] = gammap * oogammaq;
1701#endif
1702
1703 // prefactors for derivative ERI relations
1704 if (deriv_order_ > 0) {
1705#if LIBINT2_DEFINED(eri, alpha1_rho_over_zeta2)
1706 primdata.alpha1_rho_over_zeta2[0] =
1707 alpha0 * (oogammap * gammaq_o_gammapgammaq);
1708#endif
1709#if LIBINT2_DEFINED(eri, alpha2_rho_over_zeta2)
1710 primdata.alpha2_rho_over_zeta2[0] =
1711 alpha1 * (oogammap * gammaq_o_gammapgammaq);
1712#endif
1713#if LIBINT2_DEFINED(eri, alpha3_rho_over_eta2)
1714 primdata.alpha3_rho_over_eta2[0] =
1715 alpha2 * (oogammaq * gammap_o_gammapgammaq);
1716#endif
1717#if LIBINT2_DEFINED(eri, alpha4_rho_over_eta2)
1718 primdata.alpha4_rho_over_eta2[0] =
1719 alpha3 * (oogammaq * gammap_o_gammapgammaq);
1720#endif
1721#if LIBINT2_DEFINED(eri, alpha1_over_zetapluseta)
1722 primdata.alpha1_over_zetapluseta[0] = alpha0 * oogammapq;
1723#endif
1724#if LIBINT2_DEFINED(eri, alpha2_over_zetapluseta)
1725 primdata.alpha2_over_zetapluseta[0] = alpha1 * oogammapq;
1726#endif
1727#if LIBINT2_DEFINED(eri, alpha3_over_zetapluseta)
1728 primdata.alpha3_over_zetapluseta[0] = alpha2 * oogammapq;
1729#endif
1730#if LIBINT2_DEFINED(eri, alpha4_over_zetapluseta)
1731 primdata.alpha4_over_zetapluseta[0] = alpha3 * oogammapq;
1732#endif
1733#if LIBINT2_DEFINED(eri, rho12_over_alpha1)
1734 primdata.rho12_over_alpha1[0] = alpha1 * oogammap;
1735#endif
1736#if LIBINT2_DEFINED(eri, rho12_over_alpha2)
1737 primdata.rho12_over_alpha2[0] = alpha0 * oogammap;
1738#endif
1739#if LIBINT2_DEFINED(eri, rho34_over_alpha3)
1740 primdata.rho34_over_alpha3[0] = alpha3 * oogammaq;
1741#endif
1742#if LIBINT2_DEFINED(eri, rho34_over_alpha4)
1743 primdata.rho34_over_alpha4[0] = alpha2 * oogammaq;
1744#endif
1745#if LIBINT2_DEFINED(eri, two_alpha0_bra)
1746 primdata.two_alpha0_bra[0] = 2.0 * alpha0;
1747#endif
1748#if LIBINT2_DEFINED(eri, two_alpha0_ket)
1749 primdata.two_alpha0_ket[0] = 2.0 * alpha1;
1750#endif
1751#if LIBINT2_DEFINED(eri, two_alpha1_bra)
1752 primdata.two_alpha1_bra[0] = 2.0 * alpha2;
1753#endif
1754#if LIBINT2_DEFINED(eri, two_alpha1_ket)
1755 primdata.two_alpha1_ket[0] = 2.0 * alpha3;
1756#endif
1757 }
1758 } // m != 0
1759
1760 ++p;
1761 } // prefac-based prim quartet screen
1762
1763 } // rough prim quartet screen based on pair values
1764 } // ket prim pair
1765 } // bra prim pair
1766 primdata_[0].contrdepth = p;
1767 }
1768
1769#ifdef LIBINT2_ENGINE_TIMERS
1770 const auto t0 = timers.stop(0);
1771#ifdef LIBINT2_ENGINE_PROFILE_CLASS
1772 class_profiles[id].prereqs += t0.count();
1773 if (primdata_[0].contrdepth != 0) {
1774 class_profiles[id].nshellset += 1;
1775 class_profiles[id].nprimset += primdata_[0].contrdepth;
1776 }
1777#endif
1778#endif
1779
1780 // all primitive combinations screened out? set 1st target ptr to nullptr
1781 if (primdata_[0].contrdepth == 0) {
1782 targets_[0] = nullptr;
1783 return targets_;
1784 }
1785
1786 // compute directly (ss|ss)
1787 const auto compute_directly = lmax == 0 && deriv_order_ == 0;
1788
1789 if (compute_directly) {
1790#ifdef LIBINT2_ENGINE_TIMERS
1791 timers.start(1);
1792#endif
1793 auto& stack = primdata_[0].stack[0];
1794 stack = 0;
1795 for (auto p = 0; p != primdata_[0].contrdepth; ++p)
1796 stack += primdata_[p].LIBINT_T_SS_EREP_SS(0)[0];
1797 primdata_[0].targets[0] = primdata_[0].stack;
1798#ifdef LIBINT2_ENGINE_TIMERS
1799 const auto t1 = timers.stop(1);
1800#ifdef LIBINT2_ENGINE_PROFILE_CLASS
1801 class_profiles[id].build_vrr += t1.count();
1802#endif
1803#endif
1804 } // compute directly
1805 else { // call libint
1806#ifdef LIBINT2_ENGINE_TIMERS
1807#ifdef LIBINT2_PROFILE
1808 const auto t1_hrr_start = primdata_[0].timers->read(0);
1809 const auto t1_vrr_start = primdata_[0].timers->read(1);
1810#endif
1811 timers.start(1);
1812#endif
1813
1814 size_t buildfnidx;
1815 switch (braket_) {
1816 case BraKet::xx_xx:
1817 assert(bra1.contr[0].l <= hard_lmax_ &&
1818 "the angular momentum limit is exceeded");
1819 assert(bra2.contr[0].l <= hard_lmax_ &&
1820 "the angular momentum limit is exceeded");
1821 assert(ket1.contr[0].l <= hard_lmax_ &&
1822 "the angular momentum limit is exceeded");
1823 assert(ket2.contr[0].l <= hard_lmax_ &&
1824 "the angular momentum limit is exceeded");
1825 buildfnidx =
1826 ((bra1.contr[0].l * hard_lmax_ + bra2.contr[0].l) * hard_lmax_ +
1827 ket1.contr[0].l) *
1828 hard_lmax_ +
1829 ket2.contr[0].l;
1830 break;
1831
1832 case BraKet::xx_xs:
1833 assert(false && "this braket is not supported");
1834 abort();
1835 break;
1836 case BraKet::xs_xx: {
1838 int ket_lmax = hard_lmax_;
1839 switch (deriv_order_) {
1840 // N.B. notice extra PP_CAT to avoid using
1841 // LIBINT2_CENTER_DEPENDENT_MAX_AM_3eri as a subtoken which gets
1842 // expanded too early ... i.e. PP_CAT is "not associative"
1843#define BOOST_PP_NBODYENGINE_MCR8(r, data, i, elem) \
1844 case i: \
1845 BOOST_PP_IF(BOOST_PP_IS_1(BOOST_PP_CAT( \
1846 LIBINT2_CENTER_DEPENDENT_MAX_AM_, \
1847 BOOST_PP_CAT(3eri, BOOST_PP_IIF(BOOST_PP_GREATER(i, 0), i, \
1848 BOOST_PP_EMPTY())))), \
1849 ket_lmax = hard_default_lmax_, BOOST_PP_EMPTY()); \
1850 break;
1851
1852 BOOST_PP_LIST_FOR_EACH_I(BOOST_PP_NBODYENGINE_MCR8, _,
1853 BOOST_PP_NBODY_DERIV_ORDER_LIST)
1854
1855 default:
1856 assert(false && "missing case in switch");
1857 abort();
1858 }
1859 assert(bra1.contr[0].l <= hard_lmax_ &&
1860 "the angular momentum limit is exceeded");
1861 assert(ket1.contr[0].l <= ket_lmax &&
1862 "the angular momentum limit is exceeded");
1863 assert(ket2.contr[0].l <= ket_lmax &&
1864 "the angular momentum limit is exceeded");
1865 buildfnidx = (bra1.contr[0].l * ket_lmax + ket1.contr[0].l) * ket_lmax +
1866 ket2.contr[0].l;
1867#ifdef ERI3_PURE_SH
1868 if (bra1.contr[0].l > 1)
1869 assert(bra1.contr[0].pure &&
1870 "library assumes a solid harmonics shell in bra of a 3-center "
1871 "2-body int, but a cartesian shell given");
1872#endif
1873 } break;
1874
1875 case BraKet::xs_xs:
1876 assert(bra1.contr[0].l <= hard_lmax_ &&
1877 "the angular momentum limit is exceeded");
1878 assert(ket1.contr[0].l <= hard_lmax_ &&
1879 "the angular momentum limit is exceeded");
1880 buildfnidx = bra1.contr[0].l * hard_lmax_ + ket1.contr[0].l;
1881#ifdef ERI2_PURE_SH
1882 if (bra1.contr[0].l > 1)
1883 assert(bra1.contr[0].pure &&
1884 "library assumes solid harmonics shells in a 2-center "
1885 "2-body int, but a cartesian shell given in bra");
1886 if (ket1.contr[0].l > 1)
1887 assert(ket1.contr[0].pure &&
1888 "library assumes solid harmonics shells in a 2-center "
1889 "2-body int, but a cartesian shell given in bra");
1890#endif
1891 break;
1892
1893 default:
1894 assert(false && "invalid braket");
1895 abort();
1896 }
1897
1898 assert(buildfnptrs_[buildfnidx] && "null build function ptr");
1899 buildfnptrs_[buildfnidx](&primdata_[0]);
1900
1901#ifdef LIBINT2_ENGINE_TIMERS
1902 const auto t1 = timers.stop(1);
1903#ifdef LIBINT2_ENGINE_PROFILE_CLASS
1904#ifndef LIBINT2_PROFILE
1905 class_profiles[id].build_vrr += t1.count();
1906#else
1907 class_profiles[id].build_hrr += primdata_[0].timers->read(0) - t1_hrr_start;
1908 class_profiles[id].build_vrr += primdata_[0].timers->read(1) - t1_vrr_start;
1909#endif
1910#endif
1911#endif
1912
1913#ifdef LIBINT2_ENGINE_TIMERS
1914 timers.start(2);
1915#endif
1916
1917 const auto ntargets = nshellsets();
1918
1919 // if needed, permute and transform
1920 if (use_scratch) {
1921 constexpr auto using_scalar_real =
1922 sizeof(value_type) == sizeof(scalar_type);
1923 static_assert(using_scalar_real,
1924 "Libint2 C++11 API only supports scalar real types");
1925 typedef Eigen::Matrix<scalar_type, Eigen::Dynamic, Eigen::Dynamic,
1926 Eigen::RowMajor>
1927 Matrix;
1928
1929 // a 2-d view of the 4-d source tensor
1930 const auto nr1_cart = bra1.cartesian_size();
1931 const auto nr2_cart = bra2.cartesian_size();
1932 const auto nc1_cart = ket1.cartesian_size();
1933 const auto nc2_cart = ket2.cartesian_size();
1934 const auto ncol_cart = nc1_cart * nc2_cart;
1935 const auto n1234_cart = nr1_cart * nr2_cart * ncol_cart;
1936 const auto nr1 = bra1.size();
1937 const auto nr2 = bra2.size();
1938 const auto nc1 = ket1.size();
1939 const auto nc2 = ket2.size();
1940 const auto nrow = nr1 * nr2;
1941 const auto ncol = nc1 * nc2;
1942
1943 // a 2-d view of the 4-d target tensor
1944 const auto nr1_tgt = tbra1.size();
1945 const auto nr2_tgt = tbra2.size();
1946 const auto nc1_tgt = tket1.size();
1947 const auto nc2_tgt = tket2.size();
1948 const auto ncol_tgt = nc1_tgt * nc2_tgt;
1949 const auto n_tgt = nr1_tgt * nr2_tgt * ncol_tgt;
1950
1951 auto hotscr = &scratch_[0]; // points to the hot scratch
1952
1953 // transform to solid harmonics first, then unpermute, if necessary
1954 for (auto s = 0; s != ntargets; ++s) {
1955 // when permuting derivatives may need to permute shellsets also, not
1956 // just integrals
1957 // within shellsets; this will point to where source shellset s should
1958 // end up
1959 auto s_target = s;
1960
1961 auto source =
1962 primdata_[0].targets[s]; // points to the most recent result
1963 auto target = hotscr;
1964
1965 if (bra1.contr[0].pure) {
1966 libint2::solidharmonics::transform_first(
1967 bra1.contr[0].l, nr2_cart * ncol_cart, source, target);
1968 std::swap(source, target);
1969 }
1970 if (bra2.contr[0].pure) {
1971 libint2::solidharmonics::transform_inner(bra1.size(), bra2.contr[0].l,
1972 ncol_cart, source, target);
1973 std::swap(source, target);
1974 }
1975 if (ket1.contr[0].pure) {
1976 libint2::solidharmonics::transform_inner(nrow, ket1.contr[0].l,
1977 nc2_cart, source, target);
1978 std::swap(source, target);
1979 }
1980 if (ket2.contr[0].pure) {
1981 libint2::solidharmonics::transform_last(
1982 bra1.size() * bra2.size() * ket1.size(), ket2.contr[0].l, source,
1983 target);
1984 std::swap(source, target);
1985 }
1986
1987 // need to permute?
1988 if (permute) {
1989 // loop over rows of the source matrix
1990 const auto* src_row_ptr = source;
1991 auto tgt_ptr = target;
1992
1993 // if permuting derivatives ints must update their derivative index
1994 // Additional BraKet types would require adding support to
1995 // DerivMapGenerator::generate_deriv_index_map
1996 if (deriv_order_) {
1997 Tensor<size_t>& mapDerivIndex =
1998 libint2::DerivMapGenerator::instance(deriv_order_, braket_);
1999 switch (braket_) {
2000 case BraKet::xx_xx: {
2001 s_target = mapDerivIndex((size_t)swap_braket, (size_t)swap_tbra,
2002 (size_t)swap_tket, (size_t)s);
2003 } break;
2004 case BraKet::xs_xx: {
2005 assert(!swap_bra);
2006 assert(!swap_braket);
2007 s_target = mapDerivIndex((size_t)0, (size_t)0,
2008 (size_t)swap_tket, (size_t)s);
2009 } break;
2010 case BraKet::xs_xs: {
2011 assert(!swap_bra);
2012 assert(!swap_ket);
2013 assert(!swap_braket);
2014 s_target = s;
2015 } break;
2016
2017 default:
2018 assert(false &&
2019 "This braket type not yet supported for geometric "
2020 "derivatives");
2021 abort();
2022 }
2023 }
2024
2025 for (auto r1 = 0; r1 != nr1; ++r1) {
2026 for (auto r2 = 0; r2 != nr2; ++r2, src_row_ptr += ncol) {
2027 typedef Eigen::Map<const Matrix> ConstMap;
2028 typedef Eigen::Map<Matrix> Map;
2029 typedef Eigen::Map<Matrix, Eigen::Unaligned,
2030 Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>
2031 StridedMap;
2032
2033 // represent this source row as a matrix
2034 ConstMap src_blk_mat(src_row_ptr, nc1, nc2);
2035
2036 // and copy to the block of the target matrix
2037 if (swap_braket) {
2038 // if swapped bra and ket, a row of source becomes a column
2039 // of
2040 // target
2041 // source row {r1,r2} is mapped to target column {r1,r2} if
2042 // !swap_tket, else to {r2,r1}
2043 const auto tgt_col_idx =
2044 !swap_tket ? r1 * nr2 + r2 : r2 * nr1 + r1;
2045 StridedMap tgt_blk_mat(
2046 tgt_ptr + tgt_col_idx, nr1_tgt, nr2_tgt,
2047 Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>(
2048 nr2_tgt * ncol_tgt, ncol_tgt));
2049 if (swap_tbra)
2050 tgt_blk_mat = src_blk_mat.transpose();
2051 else
2052 tgt_blk_mat = src_blk_mat;
2053 } else {
2054 // source row {r1,r2} is mapped to target row {r1,r2} if
2055 // !swap_tbra, else to {r2,r1}
2056 const auto tgt_row_idx =
2057 !swap_tbra ? r1 * nr2 + r2 : r2 * nr1 + r1;
2058 Map tgt_blk_mat(tgt_ptr + tgt_row_idx * ncol, nc1_tgt, nc2_tgt);
2059 if (swap_tket)
2060 tgt_blk_mat = src_blk_mat.transpose();
2061 else
2062 tgt_blk_mat = src_blk_mat;
2063 }
2064 } // end of loop
2065 } // over rows of source
2066 std::swap(source, target);
2067 } // need to permute?
2068
2069 // if the integrals ended up in scratch_, keep them there, update the
2070 // hot buffer
2071 // to the next available scratch space, and update targets_
2072 if (source != primdata_[0].targets[s]) {
2073 hotscr += n1234_cart;
2074 if (s != s_target)
2075 assert(set_targets_ && "logic error"); // mess if targets_ points
2076 // to primdata_[0].targets
2077 targets_[s_target] = source;
2078 } else {
2079 // only needed if permuted derivs or set_targets_ is true
2080 // for simplicity always set targets_
2081 if (s != s_target)
2082 assert(set_targets_ && "logic error"); // mess if targets_ points
2083 // to primdata_[0].targets
2084 targets_[s_target] = source;
2085 }
2086 } // loop over shellsets
2087 } // if need_scratch => needed to transpose and/or tform
2088 else { // did not use scratch? may still need to update targets_
2089 if (set_targets_) {
2090 for (auto s = 0; s != ntargets; ++s)
2091 targets_[s] = primdata_[0].targets[s];
2092 }
2093 }
2094
2095#ifdef LIBINT2_ENGINE_TIMERS
2096 const auto t2 = timers.stop(2);
2097#ifdef LIBINT2_ENGINE_PROFILE_CLASS
2098 class_profiles[id].tform += t2.count();
2099#endif
2100#endif
2101 } // not (ss|ss)
2102
2103 if (cartesian_shell_normalization() == CartesianShellNormalization::uniform) {
2104 std::array<std::reference_wrapper<const Shell>, 4> shells{tbra1, tbra2,
2105 tket1, tket2};
2106 for (auto s = 0ul; s != targets_.size(); ++s) {
2107 uniform_normalize_cartesian_shells(const_cast<value_type*>(targets_[s]),
2108 shells);
2109 }
2110 }
2111
2112 return targets_;
2113}
2114
2115#undef BOOST_PP_NBODY_OPERATOR_LIST
2116#undef BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE
2117#undef BOOST_PP_NBODY_OPERATOR_INDEX_LIST
2118#undef BOOST_PP_NBODY_BRAKET_INDEX_TUPLE
2119#undef BOOST_PP_NBODY_BRAKET_INDEX_LIST
2120#undef BOOST_PP_NBODY_DERIV_ORDER_TUPLE
2121#undef BOOST_PP_NBODY_DERIV_ORDER_LIST
2122#undef BOOST_PP_NBODYENGINE_MCR3
2123#undef BOOST_PP_NBODYENGINE_MCR3_ncenter
2124#undef BOOST_PP_NBODYENGINE_MCR3_default_ncenter
2125#undef BOOST_PP_NBODYENGINE_MCR3_NCENTER
2126#undef BOOST_PP_NBODYENGINE_MCR3_OPER
2127#undef BOOST_PP_NBODYENGINE_MCR3_DERIV
2128#undef BOOST_PP_NBODYENGINE_MCR3_task
2129#undef BOOST_PP_NBODYENGINE_MCR3_TASK
2130#undef BOOST_PP_NBODYENGINE_MCR4
2131#undef BOOST_PP_NBODYENGINE_MCR5
2132#undef BOOST_PP_NBODYENGINE_MCR6
2133#undef BOOST_PP_NBODYENGINE_MCR7
2134
2135#ifdef LIBINT2_DOES_NOT_INLINE_ENGINE
2136template any Engine::enforce_params_type<Engine::empty_pod>(
2137 Operator oper, const Engine::empty_pod& params, bool throw_if_wrong_type);
2138
2139template const Engine::target_ptr_vec& Engine::compute<Shell>(
2140 const Shell& first_shell, const Shell&);
2141
2142template const Engine::target_ptr_vec& Engine::compute<Shell, Shell>(
2143 const Shell& first_shell, const Shell&, const Shell&);
2144
2145template const Engine::target_ptr_vec& Engine::compute<Shell, Shell, Shell>(
2146 const Shell& first_shell, const Shell&, const Shell&, const Shell&);
2147#endif
2148
2149} // namespace libint2
2150
2151#endif /* _libint2_src_lib_libint_engineimpl_h_ */
libint2::any
a partial C++17 std::any implementation (and less efficient than can be)
Definition any.h:67
libint2
Defaults definitions for various parameters assumed by Libint.
Definition algebra.cc:24
libint2::ScreeningMethod::Original
@ Original
standard screening method:
libint2::ScreeningMethod::Conservative
@ Conservative
conservative screening:
libint2::initialized
bool initialized()
checks if the libint has been initialized.
Definition initialize.h:68
libint2::default_params
__libint2_engine_inline libint2::any default_params(const Operator &oper)
the runtime version of operator_traits<oper>::default_params()
Definition engine.impl.h:116
libint2::uniform_normalize_cartesian_shells
void uniform_normalize_cartesian_shells(Real *intset, std::array< std::reference_wrapper< const Shell >, N > shells)
rescales cartesian Gaussians to convert from standard to uniform-normalized convention
Definition cartesian.h:117
libint2::FixedOrderedIntegerPartitionIterator
Iterates over all partitions of a non-negative integer into nonnegative integers in reverse lexicog...
Definition intpart_iter.h:74
libint2::Shell
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition shell.h:720
libint2::Shell::contr
svector< Contraction > contr
contractions
Definition shell.h:742
libint2::Shell::unit
static const Shell & unit()
Definition shell.h:906
Generated by doxygen 1.10.0