LIBINT  2.1.0-stable
engine.h
1 /*
2  * This file is a part of Libint.
3  * Copyright (C) 2004-2014 Edward F. Valeev
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU Library General Public License, version 2,
7  * as published by the Free Software Foundation.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU Library General Public License
15  * along with this program. If not, see http://www.gnu.org/licenses/.
16  *
17  */
18 
19 #ifndef _libint2_src_lib_libint_engine_h_
20 #define _libint2_src_lib_libint_engine_h_
21 
22 #include <libint2/cxxstd.h>
23 #if LIBINT2_CPLUSPLUS_STD < 2011
24 # error "libint2/engine.h requires C++11 support"
25 #endif
26 
27 #include <iostream>
28 #include <array>
29 #include <vector>
30 #include <map>
31 #include <memory>
32 
33 #include <libint2.h>
34 #include <libint2/boys.h>
35 #include <libint2/shell.h>
36 #include <libint2/timer.h>
37 #include <libint2/solidharmonics.h>
38 #include <libint2/any.h>
39 #include <libint2/util/compressed_pair.h>
40 
41 #include <libint2/boost/preprocessor.hpp>
42 
43 // extra PP macros
44 
45 #define BOOST_PP_MAKE_TUPLE_INTERNAL(z, i, last) i BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i,last))
46 #define BOOST_PP_MAKE_TUPLE(n) ( BOOST_PP_REPEAT( n , BOOST_PP_MAKE_TUPLE_INTERNAL, BOOST_PP_DEC(n) ) )
48 
49 #include <Eigen/Core>
50 
51 // the engine will be profiled by default if library was configured with --enable-profile
52 #ifdef LIBINT2_PROFILE
53 # define LIBINT2_ENGINE_TIMERS
54 // uncomment if want to profile each integral class
55 # define LIBINT2_ENGINE_PROFILE_CLASS
56 #endif
57 // uncomment if want to profile the engine even if library was configured without --enable-profile
58 //# define LIBINT2_ENGINE_TIMERS
59 
60 #ifdef __GNUC__
61 #define DEPRECATED __attribute__((deprecated))
62 #elif defined(_MSC_VER)
63 #define DEPRECATED __declspec(deprecated)
64 #else
65 #pragma message("WARNING: You need to implement DEPRECATED for this compiler")
66 #define DEPRECATED
67 #endif
68 
69 namespace libint2 {
70 
71  constexpr size_t num_geometrical_derivatives(size_t ncenter,
72  size_t deriv_order) {
73  return (deriv_order > 0) ? (num_geometrical_derivatives(ncenter, deriv_order-1) * (3*ncenter+deriv_order-1))/deriv_order : 1;
74  }
75 
76  template <typename T, unsigned N>
77  typename std::remove_all_extents<T>::type *
78  to_ptr1(T (&a)[N]) { return reinterpret_cast<
79  typename std::remove_all_extents<T>::type *
80  >(&a); }
81 
82 #if defined(LIBINT2_SUPPORT_ONEBODY)
83 
87  class OneBodyEngine {
88 
89  private:
90  typedef struct {} empty_pod;
91 
92  public:
93 
97  overlap=0,
98  kinetic=1,
99  nuclear=2,
103  _invalid=-1
104  };
105 
107 #define BOOST_PP_ONEBODY_OPERATOR_LIST (overlap, \
108  (kinetic, \
109  (elecpot, \
110  (1emultipole, \
111  (2emultipole, \
112  (3emultipole, \
113  BOOST_PP_NIL))))))
114 
115 #define BOOST_PP_ONEBODY_OPERATOR_INDEX_TUPLE BOOST_PP_MAKE_TUPLE( BOOST_PP_LIST_SIZE(BOOST_PP_ONEBODY_OPERATOR_LIST) )
116 #define BOOST_PP_ONEBODY_OPERATOR_INDEX_LIST BOOST_PP_TUPLE_TO_LIST( BOOST_PP_ONEBODY_OPERATOR_INDEX_TUPLE )
117 
118 // make list of derivative orders for 1-body ints
119 #define BOOST_PP_ONEBODY_DERIV_ORDER_TUPLE BOOST_PP_MAKE_TUPLE( BOOST_PP_INC(LIBINT2_DERIV_ONEBODY_ORDER) )
120 #define BOOST_PP_ONEBODY_DERIV_ORDER_LIST BOOST_PP_TUPLE_TO_LIST( BOOST_PP_ONEBODY_DERIV_ORDER_TUPLE )
121 
124  DEPRECATED typedef operator_type integral_type;
125 
128  template <operator_type O> struct operator_traits;
129 
131 
132 
134  OneBodyEngine() : type_(_invalid), primdata_(), stack_size_(0), lmax_(-1) {}
135 
137 
145  template <typename Params = empty_pod>
147  size_t max_nprim,
148  int max_l,
149  int deriv_order = 0,
150  Params params = empty_pod()) :
151  type_(type),
152  primdata_(max_nprim * max_nprim),
153  lmax_(max_l),
154  deriv_order_(deriv_order),
155  params_(enforce_params_type(type,params)),
156  fm_eval_(type == nuclear ? coulomb_core_eval_t::instance(2*max_l+deriv_order, 1e-25) : decltype(fm_eval_){})
157  {
158  initialize();
159  }
160 
162  // intel does not accept "move ctor = default"
164  type_(other.type_),
165  primdata_(std::move(other.primdata_)),
166  stack_size_(other.stack_size_),
167  lmax_(other.lmax_),
168  hard_lmax_(other.hard_lmax_),
169  deriv_order_(other.deriv_order_),
170  params_(std::move(other.params_)),
171  fm_eval_(std::move(other.fm_eval_)),
172  scratch_(std::move(other.scratch_)),
173  scratch2_(other.scratch2_),
174  buildfnptrs_(other.buildfnptrs_) {
175  }
176 
178  OneBodyEngine(const OneBodyEngine& other) :
179  type_(other.type_),
180  primdata_(other.primdata_.size()),
181  stack_size_(other.stack_size_),
182  lmax_(other.lmax_),
183  deriv_order_(other.deriv_order_),
184  params_(other.params_),
185  fm_eval_(other.fm_eval_) {
186  initialize();
187  }
188 
189  ~OneBodyEngine() {
190  finalize();
191  }
192 
195  type_ = other.type_;
196  primdata_ = std::move(other.primdata_);
197  stack_size_ = other.stack_size_;
198  lmax_ = other.lmax_;
199  hard_lmax_ = other.hard_lmax_;
200  deriv_order_ = other.deriv_order_;
201  params_ = std::move(other.params_);
202  fm_eval_ = std::move(other.fm_eval_);
203  scratch_ = std::move(other.scratch_);
204  scratch2_ = other.scratch2_;
205  buildfnptrs_ = other.buildfnptrs_;
206  return *this;
207  }
208 
211  type_ = other.type_;
212  primdata_.resize(other.primdata_.size());
213  stack_size_ = other.stack_size_;
214  lmax_ = other.lmax_;
215  deriv_order_ = other.deriv_order_;
216  params_ = other.params_;
217  fm_eval_ = other.fm_eval_;
218  initialize();
219  return *this;
220  }
221 
224  template <typename Params>
225  void set_params(const Params& params) {
226  params_ = params;
227  reset_scratch();
228  }
231  template <typename Params>
232  DEPRECATED void set_q(const Params& params) {
233  set_params(params);
234  }
235 
240  unsigned int nshellsets() const {
241  const unsigned int num_operator_geometrical_derivatives = (type_ == nuclear) ? this->nparams() : 0;
242  const auto ncenters = 2 + num_operator_geometrical_derivatives;
243  return nopers() * num_geometrical_derivatives(ncenters,deriv_order_);
244  }
245 
254  static unsigned int to_target_deriv_index(unsigned int deriv_idx,
255  unsigned int deriv_order,
256  unsigned int ncenters,
257  unsigned int omitted_center) {
258  auto ncenters_reduced = ncenters - 1;
259  auto nderiv_1d = 3 * ncenters_reduced;
260  assert(deriv_idx < num_geometrical_derivatives(ncenters_reduced,deriv_order));
261  switch (deriv_order) {
262  case 1: return deriv_idx >= omitted_center*3 ? deriv_idx+3 : deriv_idx;
263  default: assert(deriv_order<=1); // not implemented, won't be needed when Libint computes all derivatives
264  }
265  assert(false); // unreachable
266  return 0;
267  }
268 
271  const real_t* compute(const libint2::Shell& s1,
272  const libint2::Shell& s2) {
273 
274  // can only handle 1 contraction at a time
275  assert(s1.ncontr() == 1 && s2.ncontr() == 1);
276 
277  const auto l1 = s1.contr[0].l;
278  const auto l2 = s2.contr[0].l;
279 
280  // if want nuclear, make sure there is at least one nucleus .. otherwise the user likely forgot to call set_params
281  if (type_ == nuclear and nparams() == 0)
282  throw std::runtime_error("libint2::OneBodyEngine<nuclear>, but no charges found; forgot to call set_params()?");
283 
284  const auto n1 = s1.size();
285  const auto n2 = s2.size();
286  const auto n12 = n1 * n2;
287  const auto ncart1 = s1.cartesian_size();
288  const auto ncart2 = s2.cartesian_size();
289  const auto ncart12 = ncart1 * ncart2;
290  const auto nops = nopers();
291 
292  const auto tform_to_solids = s1.contr[0].pure || s2.contr[0].pure;
293 
294  // assert # of primitive pairs
295  const auto nprim1 = s1.nprim();
296  const auto nprim2 = s2.nprim();
297  const auto nprimpairs = nprim1 * nprim2;
298  assert(nprimpairs <= primdata_.size());
299 
300  auto nparam_sets = nparams();
301 
302  // how many shell sets will be returned?
303  auto num_shellsets = nshellsets();
304  // Libint computes derivatives with respect to one center fewer, will use translational invariance to recover
305  const auto geometry_independent_operator = type_ == overlap || type_ == kinetic;
306  const auto num_deriv_centers_computed = geometry_independent_operator ? 1 : 2;
307  auto num_shellsets_computed = nopers() *
308  num_geometrical_derivatives(num_deriv_centers_computed,
309  deriv_order_);
310  // size of ints block computed by Libint
311  const auto target_buf_size = num_shellsets_computed * ncart12;
312 
313  // will use scratch_ if:
314  // - Coulomb ints are computed 1 charge at a time, contributions are accumulated in scratch_ (unless la==lb==0)
315  // - derivatives on the missing center need to be reconstructed (no need to accumulate into scratch though)
316  // will only use scratch to accumulate ints when
317  const auto accumulate_ints_in_scratch = (type_ == nuclear);
318 
319  // adjust max angular momentum, if needed
320  const auto lmax = std::max(l1, l2);
321  assert (lmax <= lmax_);
322 
323  // where cartesian ints are located varies, sometimes we compute them in scratch, etc.
324  // this is the most likely location
325  auto cartesian_ints = primdata_[0].stack;
326 
327  // simple (s|s) ints will be computed directly and accumulated in the first element of stack
328  const auto compute_directly = lmax == 0 && deriv_order_ == 0 && (type_ == overlap || type_ == nuclear);
329  if (compute_directly) {
330  primdata_[0].stack[0] = 0;
331  }
332  else if (accumulate_ints_in_scratch)
333  memset(static_cast<void*>(&scratch_[0]), 0, sizeof(real_t)*target_buf_size);
334 
335  // loop over accumulation batches
336  for(auto pset=0u; pset!=nparam_sets; ++pset) {
337 
338  if (type_!=nuclear) assert(nparam_sets == 1);
339 
340  auto p12 = 0;
341  for(auto p1=0; p1!=nprim1; ++p1) {
342  for(auto p2=0; p2!=nprim2; ++p2, ++p12) {
343  compute_primdata(primdata_[p12],s1,s2,p1,p2,pset);
344  }
345  }
346  primdata_[0].contrdepth = p12;
347 
348  if (compute_directly) {
349  auto& result = cartesian_ints[0];
350  switch (type_) {
351  case overlap:
352  for(auto p12=0; p12 != primdata_[0].contrdepth; ++p12)
353  result += primdata_[p12]._0_Overlap_0_x[0]
354  * primdata_[p12]._0_Overlap_0_y[0]
355  * primdata_[p12]._0_Overlap_0_z[0];
356  break;
357  case nuclear:
358  for(auto p12=0; p12 != primdata_[0].contrdepth; ++p12)
359  result += primdata_[p12].LIBINT_T_S_ELECPOT_S(0)[0];
360  break;
361  default:
362  assert(false);
363  }
364  primdata_[0].targets[0] = cartesian_ints;
365  }
366  else {
367 
368  buildfnptrs_[s1.contr[0].l*hard_lmax_ + s2.contr[0].l](&primdata_[0]);
369  cartesian_ints = primdata_[0].targets[0];
370 
371  if (accumulate_ints_in_scratch) {
372  cartesian_ints = &scratch_[0];
373  std::transform(primdata_[0].targets[0], primdata_[0].targets[0] + target_buf_size,
374  &scratch_[0],
375  &scratch_[0], std::plus<real_t>());
376 
377  // need to reconstruct derivatives of nuclear ints for each nucleus
378  if (deriv_order_ > 0){
379  const auto nints_per_center = target_buf_size/2;
380  // first two blocks are derivatives with respect to Gaussian positions
381  // rest are derivs with respect to nuclear coordinates
382  auto dest = &scratch_[0] + (2+pset)*nints_per_center;
383  auto src = primdata_[0].targets[0];
384  for(auto i=0; i!=nints_per_center; ++i) {
385  dest[i] = -src[i];
386  }
387  src = primdata_[0].targets[0] + nints_per_center;
388  for(auto i=0; i!=nints_per_center; ++i) {
389  dest[i] -= src[i];
390  }
391  num_shellsets_computed+=3; // we just added 3 shell sets
392  } // reconstruct derivatives
393 
394  }
395  } // ltot != 0
396 
397  } // pset (accumulation batches)
398 
399  auto result = cartesian_ints; // will be adjusted as we proceed
400  if (tform_to_solids) {
401  // where do spherical ints go?
402  auto* spherical_ints = (cartesian_ints == &scratch_[0]) ? scratch2_ : &scratch_[0];
403  result = spherical_ints;
404 
405  // transform to solid harmonics, one shell set at a time:
406  // for each computed shell set ...
407  for(auto s=0ul; s!=num_shellsets_computed; ++s, cartesian_ints+=ncart12) {
408  // ... find its index in the target shell set:
409  // 1. if regular ints do nothing
410  // 2. for derivatives the target set includes derivatives w.r.t omitted centers,
411  // to be computed later (for all ints) or already computed (for nuclear);
412  // in the former case the "omitted" set of derivatives always comes at the end
413  // hence the index of the current shell set does not change (for 2-body ints
414  // the rules are different, but Libint will eliminate the need to reconstruct via
415  // translational invariance soon, so this logic will be unnecessary).
416  auto s_target = s;
417  // .. and compute the destination
418  spherical_ints = result + n12 * s_target;
419  if (s1.contr[0].pure && s2.contr[0].pure) {
420  libint2::solidharmonics::tform(l1, l2, cartesian_ints, spherical_ints);
421  }
422  else {
423  if (s1.contr[0].pure)
424  libint2::solidharmonics::tform_rows(l1, n2, cartesian_ints, spherical_ints);
425  else
426  libint2::solidharmonics::tform_cols(n1, l2, cartesian_ints, spherical_ints);
427  }
428  } // loop cartesian shell set
429 
430  } // tform to solids
431 
432  // if computing derivatives of ints of geometry-independent operators
433  // compute the omitted derivatives using translational invariance
434  if (deriv_order_ > 0 && geometry_independent_operator) {
435  assert(deriv_order_ == 1); // assuming 1st-order derivs here, arbitrary order later
436 
437  const auto nints_computed = n12*num_shellsets_computed; // target # of ints is twice this
438 
439  // make sure there is enough room left in libint stack
440  // if not, copy into scratch2_
441  if (not tform_to_solids) {
442  const auto stack_size_remaining = stack_size_ - (result-primdata_[0].stack) - nints_computed;
443  const auto copy_to_scratch2 = stack_size_remaining < nints_computed;
444  if (copy_to_scratch2) {
445  // this is tricky ... copy does not allow scratch2_ in [result, result + nints_computed)
446  // but this would only happen in scratch2_ == result, but definition of scratch2_ ensures this
447  std::copy(result, result + nints_computed, scratch2_);
448  result = scratch2_;
449  }
450  }
451 
452  const auto src = result;
453  const auto dest = result + nints_computed;
454  for(auto f=0ul; f!=nints_computed; ++f) {
455  dest[f] = -src[f];
456  }
457  } // rebuild omitted derivatives of Cartesian ints
458 
459  return result;
460  }
461 
462  inline void compute_primdata(Libint_t& primdata,
463  const Shell& s1, const Shell& s2,
464  size_t p1, size_t p2,
465  size_t oset);
466 
467  private:
468  operator_type type_;
469  std::vector<Libint_t> primdata_;
470  size_t stack_size_; // amount allocated by libint2_init_xxx in primdata_[0].stack
471  int lmax_;
472  int hard_lmax_; // max L supported by library for this operator type + 1
473  size_t deriv_order_;
474  any params_;
475  std::shared_ptr<coulomb_core_eval_t> fm_eval_; // this is for Coulomb only
476  std::vector<real_t> scratch_; // for transposes and/or transforming to solid harmonics
477  real_t* scratch2_; // &scratch_[0] points to the first block large enough to hold all target ints
478  // scratch2_ points to second such block. It could point into scratch_ or at primdata_[0].stack
479  typedef void (*buildfnptr_t)(const Libint_t*);
480  buildfnptr_t* buildfnptrs_;
481 
482  void reset_scratch() {
483  const auto ncart_max = (lmax_+1)*(lmax_+2)/2;
484  const auto target_shellset_size = nshellsets() * ncart_max * ncart_max;
485  // need to be able to hold 2 sets of target shellsets: the worst case occurs when dealing with
486  // 1-body Coulomb ints derivatives ... have 2+natom derivative sets that are stored in scratch
487  // then need to transform to solids. To avoid copying back and forth make sure that there is enough
488  // room to transform all ints and save them in correct order in single pass
489  const auto need_extra_large_scratch = stack_size_ < target_shellset_size;
490  scratch_.resize(need_extra_large_scratch ? 2*target_shellset_size : target_shellset_size);
491  scratch2_ = need_extra_large_scratch ? &scratch_[target_shellset_size] : primdata_[0].stack;
492  }
493 
494  void initialize() {
495  assert(libint2::initialized());
496  assert(deriv_order_ <= LIBINT2_DERIV_ONEBODY_ORDER);
497 
498 #define BOOST_PP_ONEBODYENGINE_MCR2(r,product) \
499  if (type_ == BOOST_PP_TUPLE_ELEM(2,0,product) && deriv_order_ == BOOST_PP_TUPLE_ELEM(2,1,product) ) {\
500  assert(lmax_ <= BOOST_PP_CAT(LIBINT2_MAX_AM_ , \
501  BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \
502  BOOST_PP_TUPLE_ELEM(2,0,product) ) \
503  ) ); \
504  stack_size_ = \
505  LIBINT2_PREFIXED_NAME( BOOST_PP_CAT( \
506  BOOST_PP_CAT(libint2_need_memory_ , \
507  BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \
508  BOOST_PP_TUPLE_ELEM(2,0,product) ) \
509  ), \
510  BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \
511  BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \
512  ) \
513  ) \
514  )(lmax_); \
515  LIBINT2_PREFIXED_NAME( BOOST_PP_CAT( \
516  BOOST_PP_CAT(libint2_init_ , \
517  BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \
518  BOOST_PP_TUPLE_ELEM(2,0,product) ) \
519  ), \
520  BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \
521  BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \
522  ) \
523  ) \
524  )(&primdata_[0], lmax_, 0); \
525  buildfnptrs_ = to_ptr1( \
526  LIBINT2_PREFIXED_NAME( BOOST_PP_CAT( \
527  BOOST_PP_CAT(libint2_build_ , \
528  BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \
529  BOOST_PP_TUPLE_ELEM(2,0,product) ) \
530  ), \
531  BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \
532  BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \
533  ) \
534  ) \
535  ) \
536  ); \
537  hard_lmax_ = BOOST_PP_CAT( \
538  LIBINT2_MAX_AM_ , \
539  BOOST_PP_CAT( \
540  BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \
541  BOOST_PP_TUPLE_ELEM(2,0,product) ), \
542  BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \
543  BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \
544  ) \
545  ) \
546  ) + 1; \
547  reset_scratch(); \
548  return; \
549  }
550 
551 BOOST_PP_LIST_FOR_EACH_PRODUCT ( BOOST_PP_ONEBODYENGINE_MCR2, 2, (BOOST_PP_ONEBODY_OPERATOR_INDEX_LIST, BOOST_PP_ONEBODY_DERIV_ORDER_LIST) )
552 
553  assert(false); // either deriv_order_ or type_ is wrong
554  } // initialize()
555 
556  void finalize() {
557  if (primdata_.size() != 0) {
558 
559 #define BOOST_PP_ONEBODYENGINE_MCR3(r,product) \
560  LIBINT2_PREFIXED_NAME( BOOST_PP_CAT( \
561  BOOST_PP_CAT(libint2_cleanup_ , \
562  BOOST_PP_LIST_AT(BOOST_PP_ONEBODY_OPERATOR_LIST, \
563  BOOST_PP_TUPLE_ELEM(2,0,product) ) \
564  ), \
565  BOOST_PP_IIF( BOOST_PP_GREATER(BOOST_PP_TUPLE_ELEM(2,1,product),0), \
566  BOOST_PP_TUPLE_ELEM(2,1,product), BOOST_PP_EMPTY() \
567  ) \
568  ) \
569  )(&primdata_[0]); \
570  return;
571 
572 BOOST_PP_LIST_FOR_EACH_PRODUCT ( BOOST_PP_ONEBODYENGINE_MCR3, 2, (BOOST_PP_ONEBODY_OPERATOR_INDEX_LIST, BOOST_PP_ONEBODY_DERIV_ORDER_LIST) )
573 
574  }
575  } // finalize()
576 
577  //-------
578  // utils
579  //-------
580  inline unsigned int nparams() const;
581  inline unsigned int nopers() const;
585  template <typename Params>
586  static any enforce_params_type(operator_type type,
587  const Params& params,
588  bool throw_if_wrong_type = not std::is_same<Params,empty_pod>::value);
589 
590  }; // struct OneBodyEngine
591 
592  template <OneBodyEngine::operator_type Op> struct OneBodyEngine::operator_traits {
593  typedef struct {} oper_params_type;
594  static oper_params_type default_params() { return oper_params_type{}; }
595  static constexpr unsigned int nopers = 1;
596  };
597 
600  typedef std::vector<std::pair<double, std::array<double, 3>>> oper_params_type;
601  static oper_params_type default_params() { return oper_params_type{}; }
602  static constexpr unsigned int nopers = 1;
603  };
604 
608  static oper_params_type default_params() { return oper_params_type{{0.0,0.0,0.0}}; }
609  static constexpr unsigned int nopers = 4;
610  };
614  static oper_params_type default_params() { return OneBodyEngine::operator_traits<OneBodyEngine::emultipole1>::default_params(); }
615  static constexpr unsigned int nopers = OneBodyEngine::operator_traits<OneBodyEngine::emultipole1>::nopers + 6;
616  };
620  static oper_params_type default_params() { return OneBodyEngine::operator_traits<OneBodyEngine::emultipole1>::default_params(); }
621  static constexpr unsigned int nopers = OneBodyEngine::operator_traits<OneBodyEngine::emultipole2>::nopers + 10;
622  };
623 
624  inline unsigned int OneBodyEngine::nparams() const {
625  switch (type_) {
626  case nuclear:
627  return params_.as<operator_traits<nuclear>::oper_params_type>().size();
628  default:
629  return 1;
630  }
631  return 1;
632  }
633  inline unsigned int OneBodyEngine::nopers() const {
634  switch (type_) {
635 #define BOOST_PP_ONEBODYENGINE_MCR4(r,data,i,elem) case i : return operator_traits< static_cast<OneBodyEngine::operator_type> ( i ) >::nopers;
636 BOOST_PP_LIST_FOR_EACH_I ( BOOST_PP_ONEBODYENGINE_MCR4, _, BOOST_PP_ONEBODY_OPERATOR_LIST)
637  default: break;
638  }
639  assert(false); // unreachable
640  return 0;
641  }
642  template <typename Params>
643  any OneBodyEngine::enforce_params_type(operator_type type,
644  const Params& params,
645  bool throw_if_wrong_type) {
646  any result;
647  switch(type) {
648 
649 #define BOOST_PP_ONEBODYENGINE_MCR5(r,data,i,elem) \
650  case i : \
651  if (std::is_same<Params,operator_traits< static_cast<operator_type> ( i ) >::oper_params_type>::value) \
652  result = params; \
653  else { \
654  if (throw_if_wrong_type) throw std::bad_cast(); \
655  result = operator_traits<static_cast<operator_type> ( i ) >::default_params(); \
656  } \
657  break;
658 
659 BOOST_PP_LIST_FOR_EACH_I ( BOOST_PP_ONEBODYENGINE_MCR5, _, BOOST_PP_ONEBODY_OPERATOR_LIST)
660 
661  default:
662  assert(false); // missed a case?
663  }
664  return result;
665  }
666 
667  inline void OneBodyEngine::compute_primdata(Libint_t& primdata,
668  const Shell& s1, const Shell& s2,
669  size_t p1, size_t p2,
670  size_t oset) {
671 
672  const auto& A = s1.O;
673  const auto& B = s2.O;
674 
675  const auto alpha1 = s1.alpha[p1];
676  const auto alpha2 = s2.alpha[p2];
677 
678  const auto c1 = s1.contr[0].coeff[p1];
679  const auto c2 = s2.contr[0].coeff[p2];
680 
681  const auto gammap = alpha1 + alpha2;
682  const auto oogammap = 1.0 / gammap;
683  const auto rhop_over_alpha1 = alpha2 * oogammap;
684  const auto rhop = alpha1 * rhop_over_alpha1;
685  const auto Px = (alpha1 * A[0] + alpha2 * B[0]) * oogammap;
686  const auto Py = (alpha1 * A[1] + alpha2 * B[1]) * oogammap;
687  const auto Pz = (alpha1 * A[2] + alpha2 * B[2]) * oogammap;
688  const auto AB_x = A[0] - B[0];
689  const auto AB_y = A[1] - B[1];
690  const auto AB_z = A[2] - B[2];
691  const auto AB2_x = AB_x * AB_x;
692  const auto AB2_y = AB_y * AB_y;
693  const auto AB2_z = AB_z * AB_z;
694  const auto AB2 = AB2_x + AB2_y + AB2_z;
695 
696  assert (LIBINT2_SHELLQUARTET_SET == LIBINT2_SHELLQUARTET_SET_STANDARD);
697 
698  // overlap and kinetic energy ints don't use HRR, hence VRR on both centers
699  // Coulomb potential do HRR on center 1 only
700 #if LIBINT2_DEFINED(eri,PA_x)
701  primdata.PA_x[0] = Px - A[0];
702 #endif
703 #if LIBINT2_DEFINED(eri,PA_y)
704  primdata.PA_y[0] = Py - A[1];
705 #endif
706 #if LIBINT2_DEFINED(eri,PA_z)
707  primdata.PA_z[0] = Pz - A[2];
708 #endif
709 
710  if (type_ != nuclear) {
711 
712 #if LIBINT2_DEFINED(eri,PB_x)
713  primdata.PB_x[0] = Px - B[0];
714 #endif
715 #if LIBINT2_DEFINED(eri,PB_y)
716  primdata.PB_y[0] = Py - B[1];
717 #endif
718 #if LIBINT2_DEFINED(eri,PB_z)
719  primdata.PB_z[0] = Pz - B[2];
720 #endif
721  }
722 
723  if (type_ == emultipole1 || type_ == emultipole2 || type_ == emultipole3) {
724  auto& O = params_.as<operator_traits<emultipole1>::oper_params_type>(); // same as emultipoleX
725 #if LIBINT2_DEFINED(eri,BO_x)
726  primdata.BO_x[0] = B[0] - O[0];
727 #endif
728 #if LIBINT2_DEFINED(eri,BO_y)
729  primdata.BO_y[0] = B[1] - O[1];
730 #endif
731 #if LIBINT2_DEFINED(eri,BO_z)
732  primdata.BO_z[0] = B[2] - O[2];
733 #endif
734  }
735 
736 #if LIBINT2_DEFINED(eri,oo2z)
737  primdata.oo2z[0] = 0.5*oogammap;
738 #endif
739 
740  if (type_ == nuclear) { // additional factor for electrostatic potential
741  auto& params = params_.as<operator_traits<nuclear>::oper_params_type>();
742  const auto& C = params[oset].second;
743 #if LIBINT2_DEFINED(eri,PC_x)
744  primdata.PC_x[0] = Px - C[0];
745 #endif
746 #if LIBINT2_DEFINED(eri,PC_y)
747  primdata.PC_y[0] = Py - C[1];
748 #endif
749 #if LIBINT2_DEFINED(eri,PC_z)
750  primdata.PC_z[0] = Pz - C[2];
751 #endif
752  // elecpot uses HRR
753 #if LIBINT2_DEFINED(eri,AB_x)
754  primdata.AB_x[0] = A[0] - B[0];
755 #endif
756 #if LIBINT2_DEFINED(eri,AB_y)
757  primdata.AB_y[0] = A[1] - B[1];
758 #endif
759 #if LIBINT2_DEFINED(eri,AB_z)
760  primdata.AB_z[0] = A[2] - B[2];
761 #endif
762 
763  }
764 
765  decltype(c1) sqrt_PI(1.77245385090551602729816748334);
766  const auto xyz_pfac = sqrt_PI * sqrt(oogammap);
767  const auto ovlp_ss_x = exp(- rhop * AB2_x) * xyz_pfac * c1 * c2;
768  const auto ovlp_ss_y = exp(- rhop * AB2_y) * xyz_pfac;
769  const auto ovlp_ss_z = exp(- rhop * AB2_z) * xyz_pfac;
770 
771  primdata._0_Overlap_0_x[0] = ovlp_ss_x;
772  primdata._0_Overlap_0_y[0] = ovlp_ss_y;
773  primdata._0_Overlap_0_z[0] = ovlp_ss_z;
774 
775  if (type_ == kinetic || (deriv_order_ > 0)) {
776 #if LIBINT2_DEFINED(eri,two_alpha0_bra)
777  primdata.two_alpha0_bra[0] = 2.0 * alpha1;
778 #endif
779 #if LIBINT2_DEFINED(eri,two_alpha0_ket)
780  primdata.two_alpha0_ket[0] = 2.0 * alpha2;
781 #endif
782  }
783 
784  if (type_ == nuclear) {
785 #if LIBINT2_DEFINED(eri,rho12_over_alpha1) || LIBINT2_DEFINED(eri,rho12_over_alpha2)
786  if (deriv_order_ > 0) {
787 #if LIBINT2_DEFINED(eri,rho12_over_alpha1)
788  primdata.rho12_over_alpha1[0] = rhop_over_alpha1;
789 #endif
790 #if LIBINT2_DEFINED(eri,rho12_over_alpha2)
791  primdata.rho12_over_alpha2[0] = alpha1 * oogammap;
792 #endif
793  }
794 #endif
795 #if LIBINT2_DEFINED(eri,PC_x) && LIBINT2_DEFINED(eri,PC_y) && LIBINT2_DEFINED(eri,PC_z)
796  const auto PC2 = primdata.PC_x[0] * primdata.PC_x[0] +
797  primdata.PC_y[0] * primdata.PC_y[0] +
798  primdata.PC_z[0] * primdata.PC_z[0];
799  const auto U = gammap * PC2;
800  const auto mmax = s1.contr[0].l + s2.contr[0].l + deriv_order_;
801  auto* fm_ptr = &(primdata.LIBINT_T_S_ELECPOT_S(0)[0]);
802  fm_eval_->eval(fm_ptr, U, mmax);
803 
804  decltype(U) two_o_sqrt_PI(1.12837916709551257389615890312);
805  const auto q = params_.as<operator_traits<nuclear>::oper_params_type>()[oset].first;
806  const auto pfac = - q * sqrt(gammap) * two_o_sqrt_PI * ovlp_ss_x * ovlp_ss_y * ovlp_ss_z;
807  const auto m_fence = mmax + 1;
808  for(auto m=0; m!=m_fence; ++m) {
809  fm_ptr[m] *= pfac;
810  }
811 #endif
812  }
813 
814  } // OneBodyEngine::compute_primdata()
815 
816 #undef BOOST_PP_ONEBODY_OPERATOR_LIST
817 #undef BOOST_PP_ONEBODY_OPERATOR_INDEX_TUPLE
818 #undef BOOST_PP_ONEBODY_OPERATOR_INDEX_LIST
819 #undef BOOST_PP_ONEBODY_DERIV_ORDER_TUPLE
820 #undef BOOST_PP_ONEBODY_DERIV_ORDER_LIST
821 #undef BOOST_PP_ONEBODYENGINE_MCR2
822 #undef BOOST_PP_ONEBODYENGINE_MCR3
823 #undef BOOST_PP_ONEBODYENGINE_MCR4
824 #undef BOOST_PP_ONEBODYENGINE_MCR5
825 
826 #endif // LIBINT2_SUPPORT_ONEBODY
827 
835  };
836 
839  typedef std::vector<std::pair<double,double>> ContractedGaussianGeminal;
840 
841  namespace detail {
842  template <int K> struct R12_K_G12_to_Kernel;
843  template <> struct R12_K_G12_to_Kernel<-1> {
845  };
846  template <> struct R12_K_G12_to_Kernel<0> {
847  static const MultiplicativeSphericalTwoBodyKernel value = cGTG;
848  };
849  template <> struct R12_K_G12_to_Kernel<2> {
851  };
852 
853  template <MultiplicativeSphericalTwoBodyKernel Kernel> struct TwoBodyEngineDispatcher;
854 
855  } // namespace detail
856 
857  template <MultiplicativeSphericalTwoBodyKernel Kernel> struct TwoBodyEngineTraits;
858  template <> struct TwoBodyEngineTraits<Coulomb> {
860  //typedef libint2::FmEval_Chebyshev7<double> core_eval_type;
861  //typedef libint2::FmEval_Taylor<double, 7> core_eval_type;
862  typedef struct {} oper_params_type;
863  };
864  template <> struct TwoBodyEngineTraits<cGTG> {
866  typedef ContractedGaussianGeminal oper_params_type;
867  };
869  typedef libint2::GaussianGmEval<real_t, -1> core_eval_type;
870  typedef ContractedGaussianGeminal oper_params_type;
871  };
872  template <> struct TwoBodyEngineTraits<DelcGTG_square> {
874  typedef ContractedGaussianGeminal oper_params_type;
875  };
876 
877  struct delta_gm_eval {
878  void operator()(double* Gm,
879  double rho,
880  double T,
881  int mmax) {
882  const static auto one_over_two_pi = 1.0 / (2.0 * M_PI);
883  const auto G0 = exp(-T) * rho * one_over_two_pi;
884  std::fill(Gm, Gm+mmax+1, G0);
885  }
886  };
887  template <> struct TwoBodyEngineTraits<Delta> {
888  typedef libint2::GenericGmEval<delta_gm_eval> core_eval_type; // core ints are too trivial to bother
889  typedef struct {} oper_params_type;
890  };
891 
892 #ifdef LIBINT2_SUPPORT_ERI
893 
899  template <MultiplicativeSphericalTwoBodyKernel Kernel>
901  public:
902  typedef typename libint2::TwoBodyEngineTraits<Kernel>::oper_params_type oper_params_type;
903 
905  TwoBodyEngine() : primdata_(), lmax_(-1), core_eval_pack_() {
906  set_precision(std::numeric_limits<real_t>::epsilon());
907  }
908 
910 
926  TwoBodyEngine(size_t max_nprim, int max_l,
927  int deriv_order = 0,
928  real_t precision = std::numeric_limits<real_t>::epsilon(),
929  const oper_params_type& oparams = oper_params_type()) :
930  primdata_(max_nprim * max_nprim * max_nprim * max_nprim),
931  spbra_(max_nprim), spket_(max_nprim),
932  lmax_(max_l), deriv_order_(deriv_order),
933  core_eval_pack_(core_eval_type::instance(4*lmax_ + deriv_order,
934  std::numeric_limits<real_t>::epsilon()),
935  core_eval_scratch_type(4*lmax_ + deriv_order)
936  )
937  {
938  set_precision(precision);
939  initialize();
940  init_core_ints_params(oparams);
941  }
942 
944  // intel does not support "move ctor = default"
946  primdata_(std::move(other.primdata_)),
947  spbra_(std::move(other.spbra_)), spket_(std::move(other.spket_)),
948  lmax_(other.lmax_), deriv_order_(other.deriv_order_),
949  precision_(other.precision_), ln_precision_(other.ln_precision_),
950  core_eval_pack_(std::move(other.core_eval_pack_)),
951  core_ints_params_(std::move(other.core_ints_params_)),
952  scratch_(std::move(other.scratch_)) {
953  }
954 
956  TwoBodyEngine(const TwoBodyEngine& other) :
957  primdata_(other.primdata_),
958  spbra_(other.spbra_), spket_(other.spket_),
959  lmax_(other.lmax_), deriv_order_(other.deriv_order_),
960  precision_(other.precision_), ln_precision_(other.ln_precision_),
961  core_eval_pack_(other.core_eval_pack_), core_ints_params_(other.core_ints_params_)
962  {
963  initialize();
964  }
965 
966  ~TwoBodyEngine() {
967  finalize();
968  }
969 
971  // intel does not support "move asgnmt = default"
973  primdata_ = std::move(other.primdata_);
974  lmax_ = other.lmax_;
975  deriv_order_ = other.deriv_order_;
976  precision_ = other.precision_;
977  ln_precision_ = other.ln_precision_;
978  core_eval_pack_ = std::move(other.core_eval_pack_);
979  core_ints_params_ = std::move(other.core_ints_params_);
980  scratch_ = std::move(other.scratch_);
981  return *this;
982  }
983 
986  {
987  primdata_ = other.primdata_;
988  lmax_ = other.lmax_;
989  deriv_order_ = other.deriv_order_;
990  precision_ = other.precision_;
991  ln_precision_ = other.ln_precision_;
992  core_eval_pack_ = other.core_eval_pack_;
993  core_ints_params_ = other.core_ints_params_;
994  initialize();
995  return *this;
996  }
997 
998  static bool skip_core_ints;
999 
1000 #ifdef LIBINT2_ENGINE_TIMERS
1001  Timers<3> timers; // timers[0] -> prereqs
1002  // timers[1] -> build (only meaningful if LIBINT2_PROFILE is not defined
1003  // timers[2] -> tform
1004 # ifdef LIBINT2_ENGINE_PROFILE_CLASS
1005  struct class_id {
1006  size_t l[4];
1007  template <typename I> class_id(I l0, I l1, I l2, I l3) {
1008  l[0] = l0;
1009  l[1] = l1;
1010  l[2] = l2;
1011  l[3] = l3;
1012  }
1013  bool operator<(const class_id& other) const {
1014  return ordinal(l) < ordinal(other.l);
1015  }
1016  static size_t ordinal(const size_t (&l)[4]) {
1017  return ((l[0]*LIBINT2_MAX_AM_ERI + l[1])*LIBINT2_MAX_AM_ERI + l[2])*LIBINT2_MAX_AM_ERI + l[3];
1018  }
1019  std::string to_string() const {
1020  std::ostringstream oss;
1021  oss << "(" << Shell::am_symbol(l[0]) << Shell::am_symbol(l[1])
1022  << "|" << Shell::am_symbol(l[2]) << Shell::am_symbol(l[3]) << ")";
1023  return oss.str();
1024  }
1025  };
1026  struct class_profile {
1027  double prereqs;
1028  double build_hrr;
1029  double build_vrr;
1030  double tform;
1031  size_t nshellset; // total number of shell sets
1032  size_t nprimset; // total number of primitive sets
1033  class_profile() { clear(); }
1034  class_profile(const class_profile& other) = default;
1035  void clear() {
1036  prereqs = build_hrr = build_vrr = tform = 0.;
1037  nprimset = nshellset = 0;
1038  }
1039  };
1040  std::map<class_id, class_profile> class_profiles;
1041 # endif
1042 #endif
1043 
1046  const real_t* compute(const libint2::Shell& tbra1,
1047  const libint2::Shell& tbra2,
1048  const libint2::Shell& tket1,
1049  const libint2::Shell& tket2) {
1050 
1051  //
1052  // i.e. bra and ket refer to chemists bra and ket
1053  //
1054 
1055  // can only handle 1 contraction at a time
1056  assert(tbra1.ncontr() == 1 && tbra2.ncontr() == 1 &&
1057  tket1.ncontr() == 1 && tket2.ncontr() == 1);
1058 
1059  // derivatives not supported for now
1060  assert(deriv_order_ == 0);
1061 
1062 #if LIBINT2_SHELLQUARTET_SET == LIBINT2_SHELLQUARTET_SET_STANDARD // standard angular momentum ordering
1063  auto swap_bra = (tbra1.contr[0].l < tbra2.contr[0].l);
1064  auto swap_ket = (tket1.contr[0].l < tket2.contr[0].l);
1065  auto swap_braket = (tbra1.contr[0].l + tbra2.contr[0].l > tket1.contr[0].l + tket2.contr[0].l);
1066 #else // orca angular momentum ordering
1067  auto swap_bra = (tbra1.contr[0].l > tbra2.contr[0].l);
1068  auto swap_ket = (tket1.contr[0].l > tket2.contr[0].l);
1069  auto swap_braket = (tbra1.contr[0].l + tbra2.contr[0].l < tket1.contr[0].l + tket2.contr[0].l);
1070 #endif
1071  const auto& bra1 = swap_braket ? (swap_ket ? tket2 : tket1) : (swap_bra ? tbra2 : tbra1);
1072  const auto& bra2 = swap_braket ? (swap_ket ? tket1 : tket2) : (swap_bra ? tbra1 : tbra2);
1073  const auto& ket1 = swap_braket ? (swap_bra ? tbra2 : tbra1) : (swap_ket ? tket2 : tket1);
1074  const auto& ket2 = swap_braket ? (swap_bra ? tbra1 : tbra2) : (swap_ket ? tket1 : tket2);
1075 
1076  const bool tform = tbra1.contr[0].pure || tbra2.contr[0].pure || tket1.contr[0].pure || tket2.contr[0].pure;
1077  const bool use_scratch = (swap_braket || swap_bra || swap_ket || tform);
1078 
1079  // assert # of primitive pairs
1080  auto nprim_bra1 = bra1.nprim();
1081  auto nprim_bra2 = bra2.nprim();
1082  auto nprim_ket1 = ket1.nprim();
1083  auto nprim_ket2 = ket2.nprim();
1084 
1085  // adjust max angular momentum, if needed
1086  auto lmax = std::max(std::max(bra1.contr[0].l, bra2.contr[0].l), std::max(ket1.contr[0].l, ket2.contr[0].l));
1087  assert (lmax <= lmax_);
1088  if (lmax == 0) // (ss|ss) ints will be accumulated in the first element of stack
1089  primdata_[0].stack[0] = 0;
1090 
1091 #ifdef LIBINT2_ENGINE_PROFILE_CLASS
1092  class_id id(bra1.contr[0].l, bra2.contr[0].l, ket1.contr[0].l, ket2.contr[0].l);
1093  if (class_profiles.find(id) == class_profiles.end()) {
1094  class_profile dummy;
1095  class_profiles[id] = dummy;
1096  }
1097 #endif
1098 
1099  // compute primitive data
1100 #ifdef LIBINT2_ENGINE_TIMERS
1101  timers.start(0);
1102 #endif
1103  {
1104  auto p = 0;
1105  // this is far less aggressive than should be, but proper analysis
1106  // involves both bra and ket *bases* and thus cannot be done on shell-set basis
1107  // probably ln_precision_/2 - 10 is enough
1108  spbra_.init(bra1, bra2, ln_precision_);
1109  spket_.init(ket1, ket2, ln_precision_);
1110  const auto npbra = spbra_.primpairs.size();
1111  const auto npket = spket_.primpairs.size();
1112  for(auto pb=0; pb!=npbra; ++pb) {
1113  for(auto pk=0; pk!=npket; ++pk) {
1114  if (spbra_.primpairs[pb].scr + spket_.primpairs[pk].scr > ln_precision_) {
1115  if (compute_primdata(primdata_[p],bra1,bra2,ket1,ket2,spbra_,pb,spket_,pk)) {
1116  ++p;
1117  }
1118  }
1119  }
1120  }
1121  primdata_[0].contrdepth = p;
1122  }
1123 
1124 #ifdef LIBINT2_ENGINE_TIMERS
1125  const auto t0 = timers.stop(0);
1126 # ifdef LIBINT2_ENGINE_PROFILE_CLASS
1127  class_profiles[id].prereqs += t0.count();
1128  if (primdata_[0].contrdepth != 0) {
1129  class_profiles[id].nshellset += 1;
1130  class_profiles[id].nprimset += primdata_[0].contrdepth;
1131  }
1132 # endif
1133 #endif
1134 
1135  // all primitive combinations screened out? return zeroes
1136  if (primdata_[0].contrdepth == 0) {
1137  const size_t n = bra1.size() * bra2.size() * ket1.size() * ket2.size();
1138  memset(primdata_[0].stack, 0, sizeof(real_t)*n);
1139  return primdata_[0].stack;
1140  }
1141 
1142  real_t* result = nullptr;
1143 
1144  if (lmax == 0) { // (ss|ss)
1145 #ifdef LIBINT2_ENGINE_TIMERS
1146  timers.start(1);
1147 #endif
1148  auto& stack = primdata_[0].stack[0];
1149  stack = 0;
1150  for(auto p=0; p != primdata_[0].contrdepth; ++p)
1151  stack += primdata_[p].LIBINT_T_SS_EREP_SS(0)[0];
1152  primdata_[0].targets[0] = primdata_[0].stack;
1153 #ifdef LIBINT2_ENGINE_TIMERS
1154  const auto t1 = timers.stop(1);
1155 # ifdef LIBINT2_ENGINE_PROFILE_CLASS
1156  class_profiles[id].build_vrr += t1.count();
1157 # endif
1158 #endif
1159 
1160  result = primdata_[0].targets[0];
1161  }
1162  else { // not (ss|ss)
1163 #ifdef LIBINT2_ENGINE_TIMERS
1164 # ifdef LIBINT2_PROFILE
1165  const auto t1_hrr_start = primdata_[0].timers->read(0);
1166  const auto t1_vrr_start = primdata_[0].timers->read(1);
1167 # endif
1168  timers.start(1);
1169 #endif
1170  LIBINT2_PREFIXED_NAME(libint2_build_eri)[bra1.contr[0].l][bra2.contr[0].l][ket1.contr[0].l][ket2.contr[0].l](&primdata_[0]);
1171 #ifdef LIBINT2_ENGINE_TIMERS
1172  const auto t1 = timers.stop(1);
1173 # ifdef LIBINT2_ENGINE_PROFILE_CLASS
1174 # ifndef LIBINT2_PROFILE
1175  class_profiles[id].build_vrr += t1.count();
1176 # else
1177  class_profiles[id].build_hrr += primdata_[0].timers->read(0) - t1_hrr_start;
1178  class_profiles[id].build_vrr += primdata_[0].timers->read(1) - t1_vrr_start;
1179 # endif
1180 # endif
1181 #endif
1182 
1183  result = primdata_[0].targets[0];
1184 
1185 #ifdef LIBINT2_ENGINE_TIMERS
1186  timers.start(2);
1187 #endif
1188 
1189  // if needed, permute and transform
1190  if (use_scratch) {
1191 
1192  constexpr auto using_scalar_real = std::is_same<double,real_t>::value || std::is_same<float,real_t>::value;
1193  static_assert(using_scalar_real, "Libint2 C++11 API only supports fundamental real types");
1194  typedef Eigen::Matrix<real_t, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > Matrix;
1195 
1196  // a 2-d view of the 4-d source tensor
1197  const auto nr1_cart = bra1.cartesian_size();
1198  const auto nr2_cart = bra2.cartesian_size();
1199  const auto nc1_cart = ket1.cartesian_size();
1200  const auto nc2_cart = ket2.cartesian_size();
1201  const auto ncol_cart = nc1_cart * nc2_cart;
1202  const auto nr1 = bra1.size();
1203  const auto nr2 = bra2.size();
1204  const auto nc1 = ket1.size();
1205  const auto nc2 = ket2.size();
1206  const auto nrow = nr1 * nr2;
1207  const auto ncol = nc1 * nc2;
1208 
1209  // a 2-d view of the 4-d target tensor
1210  const auto nr1_tgt = tbra1.size();
1211  const auto nr2_tgt = tbra2.size();
1212  const auto nc1_tgt = tket1.size();
1213  const auto nc2_tgt = tket2.size();
1214  const auto ncol_tgt = nc1_tgt * nc2_tgt;
1215 
1216  // transform to solid harmonics first, then unpermute, if necessary
1217  auto mainbuf = result;
1218  auto scratchbuf = &scratch_[0];
1219  if (bra1.contr[0].pure) {
1220  libint2::solidharmonics::transform_first(bra1.contr[0].l, nr2_cart*ncol_cart,
1221  mainbuf, scratchbuf);
1222  std::swap(mainbuf, scratchbuf);
1223  }
1224  if (bra2.contr[0].pure) {
1225  libint2::solidharmonics::transform_inner(bra1.size(), bra2.contr[0].l, ncol_cart,
1226  mainbuf, scratchbuf);
1227  std::swap(mainbuf, scratchbuf);
1228  }
1229  if (ket1.contr[0].pure) {
1230  libint2::solidharmonics::transform_inner(nrow, ket1.contr[0].l, nc2_cart,
1231  mainbuf, scratchbuf);
1232  std::swap(mainbuf, scratchbuf);
1233  }
1234  if (ket2.contr[0].pure) {
1235  libint2::solidharmonics::transform_last(bra1.size()*bra2.size()*ket1.size(), ket2.contr[0].l,
1236  mainbuf, scratchbuf);
1237  std::swap(mainbuf, scratchbuf);
1238  }
1239 
1240  // loop over rows of the source matrix
1241  const auto* src_row_ptr = mainbuf;
1242  auto tgt_ptr = scratchbuf;
1243  for(auto r1=0; r1!=nr1; ++r1) {
1244  for(auto r2=0; r2!=nr2; ++r2, src_row_ptr+=ncol) {
1245 
1246  typedef Eigen::Map<const Matrix> ConstMap;
1247  typedef Eigen::Map<Matrix> Map;
1248  typedef Eigen::Map<Matrix, Eigen::Unaligned, Eigen::Stride<Eigen::Dynamic,Eigen::Dynamic> > StridedMap;
1249 
1250  // represent this source row as a matrix
1251  ConstMap src_blk_mat(src_row_ptr, nc1, nc2);
1252 
1253  // and copy to the block of the target matrix
1254  if (swap_braket) {
1255  // if swapped bra and ket, a row of source becomes a column of target
1256  // source row {r1,r2} is mapped to target column {r1,r2} if !swap_ket, else to {r2,r1}
1257  const auto tgt_col_idx = !swap_ket ? r1 * nr2 + r2 : r2 * nr1 + r1;
1258  StridedMap tgt_blk_mat(tgt_ptr + tgt_col_idx, nr1_tgt, nr2_tgt,
1259  Eigen::Stride<Eigen::Dynamic,Eigen::Dynamic>(nr2_tgt*ncol_tgt,ncol_tgt));
1260  if (swap_bra)
1261  tgt_blk_mat = src_blk_mat.transpose();
1262  else
1263  tgt_blk_mat = src_blk_mat;
1264  }
1265  else {
1266  // source row {r1,r2} is mapped to target row {r1,r2} if !swap_bra, else to {r2,r1}
1267  const auto tgt_row_idx = !swap_bra ? r1 * nr2 + r2 : r2 * nr1 + r1;
1268  Map tgt_blk_mat(tgt_ptr + tgt_row_idx*ncol, nc1_tgt, nc2_tgt);
1269  if (swap_ket)
1270  tgt_blk_mat = src_blk_mat.transpose();
1271  else
1272  tgt_blk_mat = src_blk_mat;
1273  }
1274 
1275  } // end of loop
1276  } // over rows of source
1277 
1278  result = scratchbuf;
1279 
1280  } // if need_scratch => needed to transpose
1281 
1282 #ifdef LIBINT2_ENGINE_TIMERS
1283  const auto t2 = timers.stop(2);
1284 # ifdef LIBINT2_ENGINE_PROFILE_CLASS
1285  class_profiles[id].tform += t2.count();
1286 # endif
1287 #endif
1288 
1289  } // not (ss|ss)
1290 
1291  return result;
1292  }
1293 
1302  void set_precision(real_t prec) {
1303  if (prec <= 0.) {
1304  precision_ = 0.;
1305  ln_precision_ = std::numeric_limits<real_t>::lowest();
1306  }
1307  else {
1308  precision_ = prec;
1309  ln_precision_ = std::log(precision_);
1310  }
1311  }
1314  real_t precision() const {
1315  return precision_;
1316  }
1317 
1318  void print_timers() {
1319 #ifdef LIBINT2_ENGINE_TIMERS
1320  std::cout << "timers: prereq = " << timers.read(0);
1321 # ifndef LIBINT2_PROFILE // if libint's profiling was on, engine's build timer will include its overhead
1322  // do not report it, detailed profiling from libint will be printed below
1323  std::cout << " build = " << timers.read(1);
1324 # endif
1325  std::cout << " tform = " << timers.read(2) << std::endl;
1326 #endif
1327 #ifdef LIBINT2_PROFILE
1328  std::cout << "build timers: hrr = " << primdata_[0].timers->read(0)
1329  << " vrr = " << primdata_[0].timers->read(1) << std::endl;
1330 #endif
1331 #ifdef LIBINT2_ENGINE_TIMERS
1332 # ifdef LIBINT2_ENGINE_PROFILE_CLASS
1333  for(const auto& p: class_profiles) {
1334  printf("{\"%s\", %10.5lf, %10.5lf, %10.5lf, %10.5lf, %ld, %ld},\n",
1335  p.first.to_string().c_str(),
1336  p.second.prereqs,
1337  p.second.build_vrr,
1338  p.second.build_hrr,
1339  p.second.tform,
1340  p.second.nshellset,
1341  p.second.nprimset);
1342  }
1343 # endif
1344 #endif
1345  }
1346 
1347  private:
1348 
1349  inline bool compute_primdata(Libint_t& primdata,
1350  const Shell& sbra1,
1351  const Shell& sbra2,
1352  const Shell& sket1,
1353  const Shell& sket2,
1354  const ShellPair& spbra, size_t pbra,
1355  const ShellPair& spket, size_t pket);
1356 
1357  std::vector<Libint_t> primdata_;
1358  ShellPair spbra_, spket_;
1359  int lmax_;
1360  size_t deriv_order_;
1361  real_t precision_;
1362  real_t ln_precision_;
1363 
1364  typedef typename libint2::TwoBodyEngineTraits<Kernel>::core_eval_type core_eval_type;
1367  core_eval_scratch_type>
1369  core_eval_pack_type core_eval_pack_;
1370 
1371  typedef oper_params_type core_ints_params_type; // currently core ints params are always same type as operator params
1372  core_ints_params_type core_ints_params_;
1374  void init_core_ints_params(const oper_params_type& oparams);
1375 
1376  std::vector<real_t> scratch_; // for transposes and/or transforming to solid harmonics
1377 
1378  friend struct detail::TwoBodyEngineDispatcher<Kernel>;
1379 
1380  void initialize() {
1381  assert(libint2::initialized());
1382  assert(lmax_ <= LIBINT2_MAX_AM_ERI);
1383  assert(deriv_order_ <= LIBINT2_DERIV_ERI_ORDER);
1384 
1385  const auto ncart_max = (lmax_+1)*(lmax_+2)/2;
1386  const auto max_shellpair_size = ncart_max * ncart_max;
1387  const auto max_shellset_size = max_shellpair_size * max_shellpair_size;
1388 
1389  switch(deriv_order_) {
1390 
1391  case 0:
1392  libint2_init_eri(&primdata_[0], lmax_, 0);
1393  scratch_.resize(max_shellset_size);
1394  break;
1395  case 1:
1396 #if LIBINT2_DERIV_ERI_ORDER > 0
1397  libint2_init_eri1(&primdata_[0], lmax_, 0);
1398  scratch_.resize(9 * max_shellset_size);
1399 #endif
1400  break;
1401  case 2:
1402 #if LIBINT2_DERIV_ERI_ORDER > 1
1403  libint2_init_eri2(&primdata_[0], lmax_, 0);
1404  scratch_.resize(45 * max_shellset_size);
1405 #endif
1406  break;
1407  default: assert(deriv_order_ < 3);
1408  }
1409 
1410 #ifdef LIBINT2_ENGINE_TIMERS
1411  timers.set_now_overhead(25);
1412 #endif
1413 #ifdef LIBINT2_PROFILE
1414  primdata_[0].timers->set_now_overhead(25);
1415 #endif
1416  }
1417 
1418  void finalize() {
1419  if (primdata_.size() != 0) {
1420  switch(deriv_order_) {
1421 
1422  case 0:
1423  libint2_cleanup_eri(&primdata_[0]);
1424  break;
1425  case 1:
1426 #if LIBINT2_DERIV_ERI_ORDER > 0
1427  libint2_cleanup_eri1(&primdata_[0]);
1428 #endif
1429  break;
1430  case 2:
1431 #if LIBINT2_DERIV_ERI_ORDER > 1
1432  libint2_cleanup_eri2(&primdata_[0]);
1433 #endif
1434  break;
1435  }
1436  }
1437  }
1438 
1439  }; // struct TwoBodyEngine
1440 
1441  namespace detail {
1442  template <> struct TwoBodyEngineDispatcher<Coulomb> {
1443  static void core_eval(TwoBodyEngine<Coulomb>* engine,
1444  real_t* Fm,
1445  int mmax,
1446  real_t T,
1447  real_t) {
1448  engine->core_eval_pack_.first()->eval(Fm, T, mmax);
1449  }
1450  };
1451 
1452  template <>
1454  static void core_eval(TwoBodyEngine<cGTG_times_Coulomb>* engine,
1455  real_t* Gm,
1456  int mmax,
1457  real_t T,
1458  real_t rho) {
1459  engine->core_eval_pack_.first()->eval(Gm, rho, T, mmax, engine->core_ints_params_, &engine->core_eval_pack_.second());
1460  }
1461  };
1462  template <>
1464  static void core_eval(TwoBodyEngine<cGTG>* engine,
1465  real_t* Gm,
1466  int mmax,
1467  real_t T,
1468  real_t rho) {
1469  engine->core_eval_pack_.first()->eval(Gm, rho, T, mmax, engine->core_ints_params_);
1470  }
1471  };
1472  template <>
1474  static void core_eval(TwoBodyEngine<DelcGTG_square>* engine,
1475  real_t* Gm,
1476  int mmax,
1477  real_t T,
1478  real_t rho) {
1479  engine->core_eval_pack_.first()->eval(Gm, rho, T, mmax, engine->core_ints_params_);
1480  }
1481  };
1482  template <>
1484  static void core_eval(TwoBodyEngine<Delta>* engine,
1485  real_t* Gm,
1486  int mmax,
1487  real_t T,
1488  real_t rho) {
1489  engine->core_eval_pack_.first()->eval(Gm, rho, T, mmax);
1490  }
1491  };
1492  }
1493 
1494  template <MultiplicativeSphericalTwoBodyKernel Kernel>
1496 
1497  template <MultiplicativeSphericalTwoBodyKernel Kernel>
1498  inline bool TwoBodyEngine<Kernel>::compute_primdata(Libint_t& primdata,
1499  const Shell& sbra1,
1500  const Shell& sbra2,
1501  const Shell& sket1,
1502  const Shell& sket2,
1503  const ShellPair& spbra, size_t pbra,
1504  const ShellPair& spket, size_t pket) {
1505 
1506  const auto& A = sbra1.O;
1507  const auto& B = sbra2.O;
1508  const auto& C = sket1.O;
1509  const auto& D = sket2.O;
1510  const auto& AB = spbra.AB;
1511  const auto& CD = spket.AB;
1512 
1513  const auto& spbrapp = spbra.primpairs[pbra];
1514  const auto& spketpp = spket.primpairs[pket];
1515  const auto& pbra1 = spbrapp.p1;
1516  const auto& pbra2 = spbrapp.p2;
1517  const auto& pket1 = spketpp.p1;
1518  const auto& pket2 = spketpp.p2;
1519 
1520  const auto alpha0 = sbra1.alpha[pbra1];
1521  const auto alpha1 = sbra2.alpha[pbra2];
1522  const auto alpha2 = sket1.alpha[pket1];
1523  const auto alpha3 = sket2.alpha[pket2];
1524 
1525  const auto c0 = sbra1.contr[0].coeff[pbra1];
1526  const auto c1 = sbra2.contr[0].coeff[pbra2];
1527  const auto c2 = sket1.contr[0].coeff[pket1];
1528  const auto c3 = sket2.contr[0].coeff[pket2];
1529 
1530  const auto amtot = sbra1.contr[0].l + sket1.contr[0].l +
1531  sbra2.contr[0].l + sket2.contr[0].l;
1532 
1533  const auto gammap = alpha0 + alpha1;
1534  const auto oogammap = spbrapp.one_over_gamma;
1535  const auto rhop = alpha0 * alpha1 * oogammap;
1536 
1537  const auto gammaq = alpha2 + alpha3;
1538  const auto oogammaq = spketpp.one_over_gamma;
1539  const auto rhoq = alpha2 * alpha3 * oogammaq;
1540 
1541  const auto& P = spbrapp.P;
1542  const auto& Q = spketpp.P;
1543  const auto PQx = P[0] - Q[0];
1544  const auto PQy = P[1] - Q[1];
1545  const auto PQz = P[2] - Q[2];
1546  const auto PQ2 = PQx * PQx + PQy * PQy + PQz * PQz;
1547 
1548  const auto K12 = spbrapp.K * spketpp.K;
1549  decltype(K12) two_times_M_PI_to_25(34.986836655249725693); // (2 \pi)^{5/2}
1550  const auto gammapq = gammap + gammaq;
1551  const auto sqrt_gammapq = sqrt(gammapq);
1552  const auto oogammapq = 1.0 / (gammapq);
1553  auto pfac = two_times_M_PI_to_25 * K12 * sqrt_gammapq * oogammapq;
1554  pfac *= c0 * c1 * c2 * c3;
1555 
1556  if (std::abs(pfac) < precision_)
1557  return false;
1558 
1559  const auto rho = gammap * gammaq * oogammapq;
1560  const auto T = PQ2*rho;
1561  auto* fm_ptr = &(primdata.LIBINT_T_SS_EREP_SS(0)[0]);
1562  const auto mmax = amtot + deriv_order_;
1563 
1564  if (!skip_core_ints)
1565  detail::TwoBodyEngineDispatcher<Kernel>::core_eval(this, fm_ptr, mmax, T, rho);
1566 
1567  for(auto m=0; m!=mmax+1; ++m) {
1568  fm_ptr[m] *= pfac;
1569  }
1570 
1571  if (mmax == 0)
1572  return true;
1573 
1574 #if LIBINT2_DEFINED(eri,PA_x)
1575  primdata.PA_x[0] = P[0] - A[0];
1576 #endif
1577 #if LIBINT2_DEFINED(eri,PA_y)
1578  primdata.PA_y[0] = P[1] - A[1];
1579 #endif
1580 #if LIBINT2_DEFINED(eri,PA_z)
1581  primdata.PA_z[0] = P[2] - A[2];
1582 #endif
1583 #if LIBINT2_DEFINED(eri,PB_x)
1584  primdata.PB_x[0] = P[0] - B[0];
1585 #endif
1586 #if LIBINT2_DEFINED(eri,PB_y)
1587  primdata.PB_y[0] = P[1] - B[1];
1588 #endif
1589 #if LIBINT2_DEFINED(eri,PB_z)
1590  primdata.PB_z[0] = P[2] - B[2];
1591 #endif
1592 
1593 #if LIBINT2_DEFINED(eri,QC_x)
1594  primdata.QC_x[0] = Q[0] - C[0];
1595 #endif
1596 #if LIBINT2_DEFINED(eri,QC_y)
1597  primdata.QC_y[0] = Q[1] - C[1];
1598 #endif
1599 #if LIBINT2_DEFINED(eri,QC_z)
1600  primdata.QC_z[0] = Q[2] - C[2];
1601 #endif
1602 #if LIBINT2_DEFINED(eri,QD_x)
1603  primdata.QD_x[0] = Q[0] - D[0];
1604 #endif
1605 #if LIBINT2_DEFINED(eri,QD_y)
1606  primdata.QD_y[0] = Q[1] - D[1];
1607 #endif
1608 #if LIBINT2_DEFINED(eri,QD_z)
1609  primdata.QD_z[0] = Q[2] - D[2];
1610 #endif
1611 
1612 #if LIBINT2_DEFINED(eri,AB_x)
1613  primdata.AB_x[0] = AB[0];
1614 #endif
1615 #if LIBINT2_DEFINED(eri,AB_y)
1616  primdata.AB_y[0] = AB[1];
1617 #endif
1618 #if LIBINT2_DEFINED(eri,AB_z)
1619  primdata.AB_z[0] = AB[2];
1620 #endif
1621 #if LIBINT2_DEFINED(eri,BA_x)
1622  primdata.BA_x[0] = -AB[0];
1623 #endif
1624 #if LIBINT2_DEFINED(eri,BA_y)
1625  primdata.BA_y[0] = -AB[1];
1626 #endif
1627 #if LIBINT2_DEFINED(eri,BA_z)
1628  primdata.BA_z[0] = -AB[2];
1629 #endif
1630 
1631 #if LIBINT2_DEFINED(eri,CD_x)
1632  primdata.CD_x[0] = CD[0];
1633 #endif
1634 #if LIBINT2_DEFINED(eri,CD_y)
1635  primdata.CD_y[0] = CD[1];
1636 #endif
1637 #if LIBINT2_DEFINED(eri,CD_z)
1638  primdata.CD_z[0] = CD[2];
1639 #endif
1640 #if LIBINT2_DEFINED(eri,DC_x)
1641  primdata.DC_x[0] = -CD[0];
1642 #endif
1643 #if LIBINT2_DEFINED(eri,DC_y)
1644  primdata.DC_y[0] = -CD[1];
1645 #endif
1646 #if LIBINT2_DEFINED(eri,DC_z)
1647  primdata.DC_z[0] = -CD[2];
1648 #endif
1649 
1650  const auto gammap_o_gammapgammaq = oogammapq * gammap;
1651  const auto gammaq_o_gammapgammaq = oogammapq * gammaq;
1652 
1653  const auto Wx = (gammap_o_gammapgammaq * P[0] + gammaq_o_gammapgammaq * Q[0]);
1654  const auto Wy = (gammap_o_gammapgammaq * P[1] + gammaq_o_gammapgammaq * Q[1]);
1655  const auto Wz = (gammap_o_gammapgammaq * P[2] + gammaq_o_gammapgammaq * Q[2]);
1656 
1657 #if LIBINT2_DEFINED(eri,WP_x)
1658  primdata.WP_x[0] = Wx - P[0];
1659 #endif
1660 #if LIBINT2_DEFINED(eri,WP_y)
1661  primdata.WP_y[0] = Wy - P[1];
1662 #endif
1663 #if LIBINT2_DEFINED(eri,WP_z)
1664  primdata.WP_z[0] = Wz - P[2];
1665 #endif
1666 #if LIBINT2_DEFINED(eri,WQ_x)
1667  primdata.WQ_x[0] = Wx - Q[0];
1668 #endif
1669 #if LIBINT2_DEFINED(eri,WQ_y)
1670  primdata.WQ_y[0] = Wy - Q[1];
1671 #endif
1672 #if LIBINT2_DEFINED(eri,WQ_z)
1673  primdata.WQ_z[0] = Wz - Q[2];
1674 #endif
1675 #if LIBINT2_DEFINED(eri,oo2z)
1676  primdata.oo2z[0] = 0.5*oogammap;
1677 #endif
1678 #if LIBINT2_DEFINED(eri,oo2e)
1679  primdata.oo2e[0] = 0.5*oogammaq;
1680 #endif
1681 #if LIBINT2_DEFINED(eri,oo2ze)
1682  primdata.oo2ze[0] = 0.5*oogammapq;
1683 #endif
1684 #if LIBINT2_DEFINED(eri,roz)
1685  primdata.roz[0] = rho*oogammap;
1686 #endif
1687 #if LIBINT2_DEFINED(eri,roe)
1688  primdata.roe[0] = rho*oogammaq;
1689 #endif
1690 
1691  // using ITR?
1692 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x)
1693  primdata.TwoPRepITR_pfac0_0_0_x[0] = - (alpha1*AB[0] + alpha3*CD[0]) * oogammap;
1694 #endif
1695 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y)
1696  primdata.TwoPRepITR_pfac0_0_0_y[0] = - (alpha1*AB[1] + alpha3*CD[1]) * oogammap;
1697 #endif
1698 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z)
1699  primdata.TwoPRepITR_pfac0_0_0_z[0] = - (alpha1*AB[2] + alpha3*CD[2]) * oogammap;
1700 #endif
1701 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x)
1702  primdata.TwoPRepITR_pfac0_1_0_x[0] = - (alpha1*AB[0] + alpha3*CD[0]) * oogammaq;
1703 #endif
1704 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y)
1705  primdata.TwoPRepITR_pfac0_1_0_y[0] = - (alpha1*AB[1] + alpha3*CD[1]) * oogammaq;
1706 #endif
1707 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z)
1708  primdata.TwoPRepITR_pfac0_1_0_z[0] = - (alpha1*AB[2] + alpha3*CD[2]) * oogammaq;
1709 #endif
1710 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x)
1711  primdata.TwoPRepITR_pfac0_0_1_x[0] = (alpha0*AB[0] + alpha2*CD[0]) * oogammap;
1712 #endif
1713 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y)
1714  primdata.TwoPRepITR_pfac0_0_1_y[0] = (alpha0*AB[1] + alpha2*CD[1]) * oogammap;
1715 #endif
1716 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z)
1717  primdata.TwoPRepITR_pfac0_0_1_z[0] = (alpha0*AB[2] + alpha2*CD[2]) * oogammap;
1718 #endif
1719 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x)
1720  primdata.TwoPRepITR_pfac0_1_1_x[0] = (alpha0*AB[0] + alpha2*CD[0]) * oogammaq;
1721 #endif
1722 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y)
1723  primdata.TwoPRepITR_pfac0_1_1_y[0] = (alpha0*AB[1] + alpha2*CD[1]) * oogammaq;
1724 #endif
1725 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z)
1726  primdata.TwoPRepITR_pfac0_1_1_z[0] = (alpha0*AB[2] + alpha2*CD[2]) * oogammaq;
1727 #endif
1728 #if LIBINT2_DEFINED(eri,eoz)
1729  primdata.eoz[0] = gammaq*oogammap;
1730 #endif
1731 #if LIBINT2_DEFINED(eri,zoe)
1732  primdata.zoe[0] = gammap*oogammaq;
1733 #endif
1734 
1735  // prefactors for derivative ERI relations
1736  if (deriv_order_ > 0) {
1737 #if LIBINT2_DEFINED(eri,alpha1_rho_over_zeta2)
1738  primdata.alpha1_rho_over_zeta2[0] = alpha0 * (oogammap * gammaq_o_gammapgammaq);
1739 #endif
1740 #if LIBINT2_DEFINED(eri,alpha2_rho_over_zeta2)
1741  primdata.alpha2_rho_over_zeta2[0] = alpha1 * (oogammap * gammaq_o_gammapgammaq);
1742 #endif
1743 #if LIBINT2_DEFINED(eri,alpha3_rho_over_eta2)
1744  primdata.alpha3_rho_over_eta2[0] = alpha2 * (oogammaq * gammap_o_gammapgammaq);
1745 #endif
1746 #if LIBINT2_DEFINED(eri,alpha4_rho_over_eta2)
1747  primdata.alpha4_rho_over_eta2[0] = alpha3 * (oogammaq * gammap_o_gammapgammaq);
1748 #endif
1749 #if LIBINT2_DEFINED(eri,alpha1_over_zetapluseta)
1750  primdata.alpha1_over_zetapluseta[0] = alpha0 * oogammapq;
1751 #endif
1752 #if LIBINT2_DEFINED(eri,alpha2_over_zetapluseta)
1753  primdata.alpha2_over_zetapluseta[0] = alpha1 * oogammapq;
1754 #endif
1755 #if LIBINT2_DEFINED(eri,alpha3_over_zetapluseta)
1756  primdata.alpha3_over_zetapluseta[0] = alpha2 * oogammapq;
1757 #endif
1758 #if LIBINT2_DEFINED(eri,alpha4_over_zetapluseta)
1759  primdata.alpha4_over_zetapluseta[0] = alpha3 * oogammapq;
1760 #endif
1761 #if LIBINT2_DEFINED(eri,rho12_over_alpha1)
1762  primdata.rho12_over_alpha1[0] = alpha1 * oogammap;
1763 #endif
1764 #if LIBINT2_DEFINED(eri,rho12_over_alpha2)
1765  primdata.rho12_over_alpha2[0] = alpha0 * oogammap;
1766 #endif
1767 #if LIBINT2_DEFINED(eri,rho34_over_alpha3)
1768  primdata.rho34_over_alpha3[0] = alpha3 * oogammaq;
1769 #endif
1770 #if LIBINT2_DEFINED(eri,rho34_over_alpha4)
1771  primdata.rho34_over_alpha4[0] = alpha2 * oogammaq;
1772 #endif
1773 #if LIBINT2_DEFINED(eri,two_alpha0_bra)
1774  primdata.two_alpha0_bra[0] = 2.0 * alpha0;
1775 #endif
1776 #if LIBINT2_DEFINED(eri,two_alpha0_ket)
1777  primdata.two_alpha0_ket[0] = 2.0 * alpha1;
1778 #endif
1779 #if LIBINT2_DEFINED(eri,two_alpha1_bra)
1780  primdata.two_alpha1_bra[0] = 2.0 * alpha2;
1781 #endif
1782 #if LIBINT2_DEFINED(eri,two_alpha1_ket)
1783  primdata.two_alpha1_ket[0] = 2.0 * alpha3;
1784 #endif
1785  }
1786 
1787  return true;
1788  }
1789 
1790  template <>
1792  const oper_params_type& oparams) {
1793  // [g12,[- \Del^2, g12] = 2 (\Del g12) \cdot (\Del g12)
1794  // (\Del exp(-a r_12^2) \cdot (\Del exp(-b r_12^2) = 4 a b (r_{12}^2 exp(- (a+b) r_{12}^2) )
1795  // i.e. need to scale each coefficient by 4 a b
1796  const auto ng = oparams.size();
1797  core_ints_params_.reserve(ng*(ng+1)/2);
1798  for(size_t b=0; b<ng; ++b)
1799  for(size_t k=0; k<=b; ++k) {
1800  const auto gexp = oparams[b].first + oparams[k].first;
1801  const auto gcoeff = oparams[b].second * oparams[k].second * (b == k ? 1 : 2); // if a != b include ab and ba
1802  const auto gcoeff_rescaled = 4 * oparams[b].first * oparams[k].first * gcoeff;
1803  core_ints_params_.push_back(std::make_pair(gexp, gcoeff_rescaled));
1804  }
1805  }
1806 
1807  template <MultiplicativeSphericalTwoBodyKernel Kernel>
1809  const oper_params_type& oparams) {
1810  core_ints_params_ = oparams;
1811  }
1812 
1813 
1814 #endif // LIBINT2_SUPPORT_ERI
1815 
1816 } // namespace libint2
1817 
1818 #endif /* _libint2_src_lib_libint_engine_h_ */
std::vector< real_t > alpha
exponents
Definition: shell.h:100
ShellPair pre-computes shell-pair data, primitive pairs are screened to finite precision.
Definition: shell.h:302
DEPRECATED typedef operator_type integral_type
alias to operator_type for backward compatibility to pre-05/13/2015 code
Definition: engine.h:124
Definition: engine.h:830
TwoBodyEngine & operator=(TwoBodyEngine &&other)
move assignment
Definition: engine.h:972
void set_params(const Params &params)
resets operator parameters; this may be useful if need to compute Coulomb potential integrals over ba...
Definition: engine.h:225
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition: shell.h:78
TwoBodyEngine(size_t max_nprim, int max_l, int deriv_order=0, real_t precision=std::numeric_limits< real_t >::epsilon(), const oper_params_type &oparams=oper_params_type())
Constructs a (usable) TwoBodyEngine.
Definition: engine.h:926
std::array< double, 3 > oper_params_type
Cartesian coordinates of the origin with respect to which the dipole moment is defined.
Definition: engine.h:607
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
( cGTG) ( cGTG)
Definition: engine.h:833
Definition: stdarray.h:18
TwoBodyEngine & operator=(const TwoBodyEngine &other)
(deep) copy assignment
Definition: engine.h:985
unsigned int nshellsets() const
reports the number of shell sets that each call to compute() produces.
Definition: engine.h:240
some evaluators need thread-local scratch, but most don&#39;t
Definition: boys.h:1325
std::vector< std::pair< double, double > > ContractedGaussianGeminal
contracted Gaussian geminal = , represented as a vector of { , } pairs
Definition: engine.h:839
emulates boost::any
Definition: any.h:17
const real_t * compute(const libint2::Shell &tbra1, const libint2::Shell &tbra2, const libint2::Shell &tket1, const libint2::Shell &tket2)
computes shell set of integrals
Definition: engine.h:1046
describes operator sets given by OneBodyOperator
Definition: engine.h:128
OneBodyEngine::operator_traits< OneBodyEngine::emultipole1 >::oper_params_type oper_params_type
Cartesian coordinates of the origin with respect to which the multipole moments are defined...
Definition: engine.h:619
Coulomb potential due to point charges.
Definition: engine.h:99
overlap + (Cartesian) electric dipole moment, , where is relative to origin
Definition: engine.h:100
static unsigned int to_target_deriv_index(unsigned int deriv_idx, unsigned int deriv_order, unsigned int ncenters, unsigned int omitted_center)
Given the Cartesian geometric derivative index that refers to center set (0...n-1) with one center om...
Definition: engine.h:254
OneBodyEngine & operator=(const OneBodyEngine &other)
(deep) copy assignment
Definition: engine.h:210
std::array< real_t, 3 > O
origin
Definition: shell.h:102
contracted Gaussian geminal times Coulomb
Definition: engine.h:832
std::string to_string(const T &x)
Converts x to its string representation.
Definition: entity.h:71
Definition: compressed_pair.h:27
Definition: engine.h:1005
void set_precision(real_t prec)
this specifies target precision for computing the integrals.
Definition: engine.h:1302
dur_t stop(size_t t)
stops timer t
Definition: timer.h:63
OneBodyEngine & operator=(OneBodyEngine &&other)
move assignment is default
Definition: engine.h:194
TwoBodyEngine(TwoBodyEngine &&other)
move constructor
Definition: engine.h:945
overlap
Definition: engine.h:97
emultipole2 + (Cartesian) electric octupole moment,
Definition: engine.h:102
Definition: engine.h:834
const real_t * compute(const libint2::Shell &s1, const libint2::Shell &s2)
computes shell set of integrals
Definition: engine.h:271
TwoBodyEngine(const TwoBodyEngine &other)
(deep) copy constructor
Definition: engine.h:956
Definition: boys.h:1320
double read(size_t t) const
reads value (in seconds) of timer t , converted to double
Definition: timer.h:70
Definition: engine.h:1026
Definition: test_eri_rys.cc:46
Computes the Boys function, $ F_m (T) = ^1 u^{2m} (-T u^2) \, { d}u $, using 3-rd order Chebyshev int...
Definition: boys.h:239
Array idential to C++0X arrays.
Definition: stdarray_bits.h:14
electronic kinetic energy, i.e.
Definition: engine.h:98
OneBodyEngine(operator_type type, size_t max_nprim, int max_l, int deriv_order=0, Params params=empty_pod())
Constructs a (usable) OneBodyEngine.
Definition: engine.h:146
emultipole1 + (Cartesian) electric quadrupole moment,
Definition: engine.h:101
MultiplicativeSphericalTwoBodyKernel
types of multiplicative spherically-symmetric two-body kernels known by TwoBodyEngine ...
Definition: engine.h:829
void start(size_t t)
starts timer t
Definition: timer.h:58
operator_type
types of operators (operator sets) supported by OneBodyEngine.
Definition: engine.h:96
Definition: engine.h:842
contracted Gaussian geminal =
Definition: engine.h:831
DEPRECATED void set_q(const Params &params)
alias to set_params() for backward compatibility with pre-05/13/2015 code
Definition: engine.h:232
OneBodyEngine::operator_traits< OneBodyEngine::emultipole1 >::oper_params_type oper_params_type
Cartesian coordinates of the origin with respect to which the multipole moments are defined...
Definition: engine.h:613
TwoBodyEngine()
creates a default (unusable) TwoBodyEngine
Definition: engine.h:905
OneBodyEngine(OneBodyEngine &&other)
move constructor
Definition: engine.h:163
std::vector< std::pair< double, std::array< double, 3 > > > oper_params_type
point charges and their positions
Definition: engine.h:600
OneBodyEngine()
creates a default (unusable) OneBodyEngine; to be used as placeholder for copying a usable engine ...
Definition: engine.h:134
Definition: engine.h:877
void set_now_overhead(size_t ns)
use this to report the overhead of now() call; if set, the reported timings will be adjusted for this...
Definition: timer.h:53
TwoBodyEngine computes (ab|O|cd) (i.e.
Definition: engine.h:900
Definition: engine.h:857
Computes the Boys function, $ F_m (T) = ^1 u^{2m} (-T u^2) \, { d}u $, using Taylor interpolation of ...
Definition: boys.h:752
OneBodyEngine(const OneBodyEngine &other)
(deep) copy constructor
Definition: engine.h:178
real_t precision() const
Definition: engine.h:1314
std::vector< Contraction > contr
contractions
Definition: shell.h:101
OneBodyEngine computes integrals of operators (or operator sets) given by OneBodyOperator::operator_t...
Definition: engine.h:87