LIBINT  2.1.0-stable
prep_libint2.h
1 /*
2  * This file is a part of Libint.
3  * Copyright (C) 2004-2014 Edward F. Valeev
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 2 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see http://www.gnu.org/licenses/.
17  *
18  */
19 
20 #ifndef PREP_LIBINT2_SKIP_BOOST
21 # include <boost/random.hpp>
22 #endif
23 
24 #include <cmath>
25 #include <vector>
26 #include <cstdlib>
27 #include <libint2.h>
28 #include <libint2/boys.h>
29 #include <time.h>
30 
31 extern libint2::FmEval_Chebyshev3<double> fmeval_chebyshev;
32 extern libint2::FmEval_Taylor<double,6> fmeval_taylor;
33 
34 #ifdef PREP_LIBINT2_SKIP_BOOST
35 struct RandomDie {
36  RandomDie(double lower, double upper, long seed = -1) : lower_(lower), frange_(upper - lower), seed_(seed)
37  {
38  if (seed == -1)
39  srandom(time(NULL));
40  else
41  srandom((unsigned int)seed);
42  }
43  double operator()() const {
44  return lower_ + random() * frange_ / intrange_;
45  }
46  double lower_, frange_;
47  unsigned int seed_;
48  static const long intrange_ = (1l<<31) - 1;
49 };
50 #endif
51 
52 typedef unsigned int uint;
53 template <unsigned int N>
55  public:
56  RandomShellSet(uint* am,
57  uint veclen,
58  uint contrdepth) {
59 
60  std::copy(am, am+N, l);
61 
62 #ifndef PREP_LIBINT2_SKIP_BOOST
63  boost::mt19937 rng; // produces randomness out of thin air
64  // rng.seed(1352054452); // uncomment to make results repeatable
65  boost::uniform_real<> rdist(0.1, 3.0); // distribution that maps to 0.1 .. 3.0
66  boost::variate_generator<boost::mt19937&, boost::uniform_real<> >
67  die(rng, rdist); // glues randomness with mapping
68 #else
69  RandomDie die(0.1, 3.0);
70 #endif
71 
72  for(uint c=0; c<N; ++c) {
73 
74  R[c].resize(3); generate(R[c].begin(), R[c].end(), die);
75 
76  exp[c].resize(veclen);
77  coef[c].resize(veclen);
78  for(uint v=0; v<veclen; ++v) {
79  exp[c][v].resize(contrdepth); generate(exp[c][v].begin(), exp[c][v].end(), die);
80  coef[c][v].resize(contrdepth); generate(coef[c][v].begin(), coef[c][v].end(), die);
81  }
82  }
83 
84  }
85 
86  uint l[N]; // angular momenta
87  std::vector<double> R[N]; // origins
88  std::vector< std::vector<double> > exp[N]; // exponents
89  std::vector< std::vector<double> > coef[N]; // coefficients
90 };
91 
92 template<typename LibintEval>
93 void prep_libint2(LibintEval* erievals,
94  const RandomShellSet<4u>& rsqset,
95  int norm_flag,
96  int deriv_order = 0) {
97  const uint veclen = rsqset.exp[0].size();
98  const uint contrdepth = rsqset.exp[0][0].size();
99 
100  const double* A = &(rsqset.R[0][0]);
101  const double* B = &(rsqset.R[1][0]);
102  const double* C = &(rsqset.R[2][0]);
103  const double* D = &(rsqset.R[3][0]);
104 
105  const uint* am = rsqset.l;
106  const unsigned int am01 = am[0] + am[1];
107  const unsigned int am23 = am[2] + am[3];
108  const unsigned int amtot = am01 + am23 + deriv_order;
109  // static seems to be important for performance on OS X 10.7 with Apple clang 3.1 (318.0.58)
110  // 8/24/2012 also seems to affect other compilers (macports g++ 4.7.1) on OS X 10.7 and 10.8
111  static double F[LIBINT_MAX_AM*4 + 6];
112 
113  uint p0123 = 0;
114  for (uint p0 = 0; p0 < contrdepth; p0++) {
115  for (uint p1 = 0; p1 < contrdepth; p1++) {
116  for (uint p2 = 0; p2 < contrdepth; p2++) {
117  for (uint p3 = 0; p3 < contrdepth; p3++, p0123++) {
118 
119  LibintEval* erieval = &erievals[p0123];
120  erieval->veclen = veclen;
121 #if LIBINT2_FLOP_COUNT
122  erieval->nflops = erievals[0].nflops;
123 #endif
124 
125  for (uint v = 0; v < veclen; v++) {
126 
127  const double alpha0 = rsqset.exp[0][v][p0];
128  const double alpha1 = rsqset.exp[1][v][p1];
129  const double alpha2 = rsqset.exp[2][v][p2];
130  const double alpha3 = rsqset.exp[3][v][p3];
131 
132  const double c0 = rsqset.coef[0][v][p0];
133  const double c1 = rsqset.coef[1][v][p1];
134  const double c2 = rsqset.coef[2][v][p2];
135  const double c3 = rsqset.coef[3][v][p3];
136 
137  const double gammap = alpha0 + alpha1;
138  const double oogammap = 1.0 / gammap;
139  const double rhop = alpha0 * alpha1 * oogammap;
140  const double Px = (alpha0 * A[0] + alpha1 * B[0]) * oogammap;
141  const double Py = (alpha0 * A[1] + alpha1 * B[1]) * oogammap;
142  const double Pz = (alpha0 * A[2] + alpha1 * B[2]) * oogammap;
143  const double AB_x = A[0] - B[0];
144  const double AB_y = A[1] - B[1];
145  const double AB_z = A[2] - B[2];
146  const double AB2 = AB_x * AB_x + AB_y * AB_y + AB_z * AB_z;
147 
148  const double PAx = Px - A[0];
149  const double PAy = Py - A[1];
150  const double PAz = Pz - A[2];
151  const double PBx = Px - B[0];
152  const double PBy = Py - B[1];
153  const double PBz = Pz - B[2];
154 
155 #if LIBINT2_DEFINED(eri,PA_x)
156  erieval->PA_x[v] = PAx;
157 #endif
158 #if LIBINT2_DEFINED(eri,PA_y)
159  erieval->PA_y[v] = PAy;
160 #endif
161 #if LIBINT2_DEFINED(eri,PA_z)
162  erieval->PA_z[v] = PAz;
163 #endif
164 #if LIBINT2_DEFINED(eri,PB_x)
165  erieval->PB_x[v] = PBx;
166 #endif
167 #if LIBINT2_DEFINED(eri,PB_y)
168  erieval->PB_y[v] = PBy;
169 #endif
170 #if LIBINT2_DEFINED(eri,PB_z)
171  erieval->PB_z[v] = PBz;
172 #endif
173 
174 #if LIBINT2_DEFINED(eri,AB_x)
175  erieval->AB_x[v] = AB_x;
176 #endif
177 #if LIBINT2_DEFINED(eri,AB_y)
178  erieval->AB_y[v] = AB_y;
179 #endif
180 #if LIBINT2_DEFINED(eri,AB_z)
181  erieval->AB_z[v] = AB_z;
182 #endif
183 #if LIBINT2_DEFINED(eri,BA_x)
184  erieval->BA_x[v] = -AB_x;
185 #endif
186 #if LIBINT2_DEFINED(eri,BA_y)
187  erieval->BA_y[v] = -AB_y;
188 #endif
189 #if LIBINT2_DEFINED(eri,BA_z)
190  erieval->BA_z[v] = -AB_z;
191 #endif
192 #if LIBINT2_DEFINED(eri,oo2z)
193  erieval->oo2z[v] = 0.5*oogammap;
194 #endif
195 
196  const double gammaq = alpha2 + alpha3;
197  const double oogammaq = 1.0 / gammaq;
198  const double rhoq = alpha2 * alpha3 * oogammaq;
199  const double one_o_gammap_plus_gammaq = 1.0 / (gammap + gammaq);
200  const double gammapq = gammap * gammaq * one_o_gammap_plus_gammaq;
201  const double gammap_o_gammapgammaq = gammapq * oogammaq;
202  const double gammaq_o_gammapgammaq = gammapq * oogammap;
203  const double Qx = (alpha2 * C[0] + alpha3 * D[0]) * oogammaq;
204  const double Qy = (alpha2 * C[1] + alpha3 * D[1]) * oogammaq;
205  const double Qz = (alpha2 * C[2] + alpha3 * D[2]) * oogammaq;
206  const double CD_x = C[0] - D[0];
207  const double CD_y = C[1] - D[1];
208  const double CD_z = C[2] - D[2];
209  const double CD2 = CD_x * CD_x + CD_y * CD_y + CD_z * CD_z;
210 
211  const double PQx = Px - Qx;
212  const double PQy = Py - Qy;
213  const double PQz = Pz - Qz;
214  const double PQ2 = PQx * PQx + PQy * PQy + PQz * PQz;
215 
216  const double QCx = Qx - C[0];
217  const double QCy = Qy - C[1];
218  const double QCz = Qz - C[2];
219  const double QDx = Qx - D[0];
220  const double QDy = Qy - D[1];
221  const double QDz = Qz - D[2];
222 
223 #if LIBINT2_DEFINED(eri,QC_x)
224  erieval->QC_x[v] = QCx;
225 #endif
226 #if LIBINT2_DEFINED(eri,QC_y)
227  erieval->QC_y[v] = QCy;
228 #endif
229 #if LIBINT2_DEFINED(eri,QC_z)
230  erieval->QC_z[v] = QCz;
231 #endif
232 #if LIBINT2_DEFINED(eri,QD_x)
233  erieval->QD_x[v] = QDx;
234 #endif
235 #if LIBINT2_DEFINED(eri,QD_y)
236  erieval->QD_y[v] = QDy;
237 #endif
238 #if LIBINT2_DEFINED(eri,QD_z)
239  erieval->QD_z[v] = QDz;
240 #endif
241 
242 #if LIBINT2_DEFINED(eri,CD_x)
243  erieval->CD_x[v] = CD_x;
244 #endif
245 #if LIBINT2_DEFINED(eri,CD_y)
246  erieval->CD_y[v] = CD_y;
247 #endif
248 #if LIBINT2_DEFINED(eri,CD_z)
249  erieval->CD_z[v] = CD_z;
250 #endif
251 #if LIBINT2_DEFINED(eri,DC_x)
252  erieval->DC_x[v] = -CD_x;
253 #endif
254 #if LIBINT2_DEFINED(eri,DC_y)
255  erieval->DC_y[v] = -CD_y;
256 #endif
257 #if LIBINT2_DEFINED(eri,DC_z)
258  erieval->DC_z[v] = -CD_z;
259 #endif
260 #if LIBINT2_DEFINED(eri,oo2e)
261  erieval->oo2e[v] = 0.5*oogammaq;
262 #endif
263 
264  // Prefactors for interelectron transfer relation
265 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x)
266  erieval->TwoPRepITR_pfac0_0_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammap;
267 #endif
268 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y)
269  erieval->TwoPRepITR_pfac0_0_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammap;
270 #endif
271 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z)
272  erieval->TwoPRepITR_pfac0_0_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammap;
273 #endif
274 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x)
275  erieval->TwoPRepITR_pfac0_1_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammaq;
276 #endif
277 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y)
278  erieval->TwoPRepITR_pfac0_1_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammaq;
279 #endif
280 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z)
281  erieval->TwoPRepITR_pfac0_1_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammaq;
282 #endif
283 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x)
284  erieval->TwoPRepITR_pfac0_0_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammap;
285 #endif
286 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y)
287  erieval->TwoPRepITR_pfac0_0_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammap;
288 #endif
289 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z)
290  erieval->TwoPRepITR_pfac0_0_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammap;
291 #endif
292 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x)
293  erieval->TwoPRepITR_pfac0_1_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammaq;
294 #endif
295 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y)
296  erieval->TwoPRepITR_pfac0_1_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammaq;
297 #endif
298 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z)
299  erieval->TwoPRepITR_pfac0_1_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammaq;
300 #endif
301 #if LIBINT2_DEFINED(eri,eoz)
302  erieval->eoz[v] = gammaq*oogammap;
303 #endif
304 #if LIBINT2_DEFINED(eri,zoe)
305  erieval->zoe[v] = gammap*oogammaq;
306 #endif
307 
308  const double Wx = (gammap_o_gammapgammaq * Px + gammaq_o_gammapgammaq * Qx);
309  const double Wy = (gammap_o_gammapgammaq * Py + gammaq_o_gammapgammaq * Qy);
310  const double Wz = (gammap_o_gammapgammaq * Pz + gammaq_o_gammapgammaq * Qz);
311 
312 #if LIBINT2_DEFINED(eri,WP_x)
313  erieval->WP_x[v] = Wx - Px;
314 #endif
315 #if LIBINT2_DEFINED(eri,WP_y)
316  erieval->WP_y[v] = Wy - Py;
317 #endif
318 #if LIBINT2_DEFINED(eri,WP_z)
319  erieval->WP_z[v] = Wz - Pz;
320 #endif
321 #if LIBINT2_DEFINED(eri,WQ_x)
322  erieval->WQ_x[v] = Wx - Qx;
323 #endif
324 #if LIBINT2_DEFINED(eri,WQ_y)
325  erieval->WQ_y[v] = Wy - Qy;
326 #endif
327 #if LIBINT2_DEFINED(eri,WQ_z)
328  erieval->WQ_z[v] = Wz - Qz;
329 #endif
330 #if LIBINT2_DEFINED(eri,oo2ze)
331  erieval->oo2ze[v] = 0.5/(gammap+gammaq);
332 #endif
333 #if LIBINT2_DEFINED(eri,roz)
334  erieval->roz[v] = gammapq*oogammap;
335 #endif
336 #if LIBINT2_DEFINED(eri,roe)
337  erieval->roe[v] = gammapq*oogammaq;
338 #endif
339 
340  if (deriv_order > 0) {
341  // prefactors for derivative ERI relations
342 #if LIBINT2_DEFINED(eri,alpha1_rho_over_zeta2)
343  erieval->alpha1_rho_over_zeta2[v] = alpha0 * gammapq / (gammap * gammap);
344 #endif
345 #if LIBINT2_DEFINED(eri,alpha2_rho_over_zeta2)
346  erieval->alpha2_rho_over_zeta2[v] = alpha1 * gammapq / (gammap * gammap);
347 #endif
348 #if LIBINT2_DEFINED(eri,alpha3_rho_over_eta2)
349  erieval->alpha3_rho_over_eta2[v] = alpha2 * gammapq / (gammaq * gammaq);
350 #endif
351 #if LIBINT2_DEFINED(eri,alpha4_rho_over_eta2)
352  erieval->alpha4_rho_over_eta2[v] = alpha3 * gammapq / (gammaq * gammaq);
353 #endif
354 #if LIBINT2_DEFINED(eri,alpha1_over_zetapluseta)
355  erieval->alpha1_over_zetapluseta[v] = alpha0 / (gammap + gammaq);
356 #endif
357 #if LIBINT2_DEFINED(eri,alpha2_over_zetapluseta)
358  erieval->alpha2_over_zetapluseta[v] = alpha1 / (gammap + gammaq);
359 #endif
360 #if LIBINT2_DEFINED(eri,alpha3_over_zetapluseta)
361  erieval->alpha3_over_zetapluseta[v] = alpha2 / (gammap + gammaq);
362 #endif
363 #if LIBINT2_DEFINED(eri,alpha4_over_zetapluseta)
364  erieval->alpha4_over_zetapluseta[v] = alpha3 / (gammap + gammaq);
365 #endif
366 #if LIBINT2_DEFINED(eri,rho12_over_alpha1)
367  erieval->rho12_over_alpha1[v] = rhop / alpha0;
368 #endif
369 #if LIBINT2_DEFINED(eri,rho12_over_alpha2)
370  erieval->rho12_over_alpha2[v] = rhop / alpha1;
371 #endif
372 #if LIBINT2_DEFINED(eri,rho34_over_alpha3)
373  erieval->rho34_over_alpha3[v] = rhoq / alpha2;
374 #endif
375 #if LIBINT2_DEFINED(eri,rho34_over_alpha4)
376  erieval->rho34_over_alpha4[v] = rhoq / alpha3;
377 #endif
378 #if LIBINT2_DEFINED(eri,two_alpha0_bra)
379  erieval->two_alpha0_bra[v] = 2.0 * alpha0;
380 #endif
381 #if LIBINT2_DEFINED(eri,two_alpha0_ket)
382  erieval->two_alpha0_ket[v] = 2.0 * alpha1;
383 #endif
384 #if LIBINT2_DEFINED(eri,two_alpha1_bra)
385  erieval->two_alpha1_bra[v] = 2.0 * alpha2;
386 #endif
387 #if LIBINT2_DEFINED(eri,two_alpha1_ket)
388  erieval->two_alpha1_ket[v] = 2.0 * alpha3;
389 #endif
390  }
391 
392  const double K1 = exp(- rhop * AB2) * oogammap;
393  const double K2 = exp(- rhoq * CD2) * oogammaq;
394  const double two_times_M_PI_to_25 = 34.986836655249725693;
395  double pfac = two_times_M_PI_to_25 * K1 * K2 * sqrt(one_o_gammap_plus_gammaq);
396  pfac *= c0 * c1 * c2 * c3;
397 
398  // if veclen=1, ignore the F scratch, use the erieval directly
399  if (veclen == 1 && std::is_same<double,LIBINT2_REALTYPE>::value) {
400  { // this is only used for double realtype, so nothing nefarious here, just a workaround the type system
401  double* ssss_ptr = reinterpret_cast<double*>(erieval->LIBINT_T_SS_EREP_SS(0));
402  fmeval_chebyshev.eval(ssss_ptr,PQ2*gammapq,amtot);
403  }
404  LIBINT2_REALTYPE* ssss_ptr = erieval->LIBINT_T_SS_EREP_SS(0);
405  for(unsigned int l=0; l<=amtot; ++l, ++ssss_ptr)
406  *ssss_ptr = *ssss_ptr * pfac;
407  }
408  else {
409  fmeval_chebyshev.eval(F,PQ2*gammapq,amtot);
410  //fmeval_taylor.eval(F,PQ2*gammapq,amtot);
411  {
412  LIBINT2_REALTYPE* ssss_ptr = erieval->LIBINT_T_SS_EREP_SS(0) + v;
413  for(unsigned int l=0; l<=amtot; ++l, ssss_ptr+=veclen)
414  *ssss_ptr = pfac*F[l];
415  }
416  }
417 
418  } // end of v loop
419 
420  }
421  }
422  }
423  } // end of primitive loops
424 }
425 
426 template<typename LibintEval>
427 void prep_libint2(LibintEval* erievals,
428  const RandomShellSet<3u>& rsqset,
429  int norm_flag,
430  int deriv_order = 0) {
431  const uint veclen = rsqset.exp[0].size();
432  const uint contrdepth = rsqset.exp[0][0].size();
433 
434  const unsigned int dummy_center = (LIBINT_SHELL_SET == LIBINT_SHELL_SET_STANDARD) ? 1 : 0;
435  const double* A = &(rsqset.R[0][0]);
436  const double* B = &(rsqset.R[0][0]);
437  const double* C = &(rsqset.R[1][0]);
438  const double* D = &(rsqset.R[2][0]);
439 
440  const uint* am = rsqset.l;
441  const unsigned int amtot = am[0] + am[1] + am[2] + deriv_order;
442  static double F[LIBINT_MAX_AM*3 + 6];
443 
444  uint p012 = 0;
445  for (uint p0 = 0; p0 < contrdepth; p0++) {
446  for (uint p1 = 0; p1 < contrdepth; p1++) {
447  for (uint p2 = 0; p2 < contrdepth; p2++, p012++) {
448 
449  LibintEval* erieval = &erievals[p012];
450  erieval->veclen = veclen;
451 #if LIBINT2_FLOP_COUNT
452  erieval->nflops = erievals[0].nflops;
453 #endif
454 
455  for (uint v = 0; v < veclen; v++) {
456 
457  const double alpha0 = (dummy_center == 0) ? 0.0 : rsqset.exp[0][v][p0];
458  const double alpha1 = (dummy_center == 1) ? 0.0 : rsqset.exp[0][v][p0];
459  const double alpha2 = rsqset.exp[1][v][p1];
460  const double alpha3 = rsqset.exp[2][v][p2];
461 
462  const double c0 = (dummy_center == 0) ? 1.0 : rsqset.coef[0][v][p0];
463  const double c1 = (dummy_center == 1) ? 1.0 : rsqset.coef[0][v][p0];
464  const double c2 = rsqset.coef[1][v][p1];
465  const double c3 = rsqset.coef[2][v][p2];
466 
467  const double gammap = alpha0 + alpha1;
468  const double oogammap = 1.0 / gammap;
469  const double rhop = alpha0 * alpha1 * oogammap;
470  const double Px = (alpha0 * A[0] + alpha1 * B[0]) * oogammap;
471  const double Py = (alpha0 * A[1] + alpha1 * B[1]) * oogammap;
472  const double Pz = (alpha0 * A[2] + alpha1 * B[2]) * oogammap;
473  const double PAx = Px - A[0];
474  const double PAy = Py - A[1];
475  const double PAz = Pz - A[2];
476  const double PBx = Px - B[0];
477  const double PBy = Py - B[1];
478  const double PBz = Pz - B[2];
479  const double AB_x = A[0] - B[0];
480  const double AB_y = A[1] - B[1];
481  const double AB_z = A[2] - B[2];
482  const double AB2 = AB_x * AB_x + AB_y * AB_y + AB_z * AB_z;
483 
484 #if LIBINT2_DEFINED(eri,PA_x)
485  erieval->PA_x[v] = PAx;
486 #endif
487 #if LIBINT2_DEFINED(eri,PA_y)
488  erieval->PA_y[v] = PAy;
489 #endif
490 #if LIBINT2_DEFINED(eri,PA_z)
491  erieval->PA_z[v] = PAz;
492 #endif
493 #if LIBINT2_DEFINED(eri,PB_x)
494  erieval->PB_x[v] = PBx;
495 #endif
496 #if LIBINT2_DEFINED(eri,PB_y)
497  erieval->PB_y[v] = PBy;
498 #endif
499 #if LIBINT2_DEFINED(eri,PB_z)
500  erieval->PB_z[v] = PBz;
501 #endif
502 
503 #if LIBINT2_DEFINED(eri,AB_x)
504  erieval->AB_x[v] = AB_x;
505 #endif
506 #if LIBINT2_DEFINED(eri,AB_y)
507  erieval->AB_y[v] = AB_y;
508 #endif
509 #if LIBINT2_DEFINED(eri,AB_z)
510  erieval->AB_z[v] = AB_z;
511 #endif
512 #if LIBINT2_DEFINED(eri,BA_x)
513  erieval->BA_x[v] = -AB_x;
514 #endif
515 #if LIBINT2_DEFINED(eri,BA_y)
516  erieval->BA_y[v] = -AB_y;
517 #endif
518 #if LIBINT2_DEFINED(eri,BA_z)
519  erieval->BA_z[v] = -AB_z;
520 #endif
521 #if LIBINT2_DEFINED(eri,oo2z)
522  erieval->oo2z[v] = 0.5*oogammap;
523 #endif
524 
525  const double gammaq = alpha2 + alpha3;
526  const double oogammaq = 1.0 / gammaq;
527  const double rhoq = alpha2 * alpha3 * oogammaq;
528  const double gammapq = gammap * gammaq / (gammap + gammaq);
529  const double gammap_o_gammapgammaq = gammapq * oogammaq;
530  const double gammaq_o_gammapgammaq = gammapq * oogammap;
531  const double Qx = (alpha2 * C[0] + alpha3 * D[0]) * oogammaq;
532  const double Qy = (alpha2 * C[1] + alpha3 * D[1]) * oogammaq;
533  const double Qz = (alpha2 * C[2] + alpha3 * D[2]) * oogammaq;
534  const double QCx = Qx - C[0];
535  const double QCy = Qy - C[1];
536  const double QCz = Qz - C[2];
537  const double QDx = Qx - D[0];
538  const double QDy = Qy - D[1];
539  const double QDz = Qz - D[2];
540  const double CD_x = C[0] - D[0];
541  const double CD_y = C[1] - D[1];
542  const double CD_z = C[2] - D[2];
543  const double CD2 = CD_x * CD_x + CD_y * CD_y + CD_z * CD_z;
544 
545 #if LIBINT2_DEFINED(eri,QC_x)
546  erieval->QC_x[v] = QCx;
547 #endif
548 #if LIBINT2_DEFINED(eri,QC_y)
549  erieval->QC_y[v] = QCy;
550 #endif
551 #if LIBINT2_DEFINED(eri,QC_z)
552  erieval->QC_z[v] = QCz;
553 #endif
554 #if LIBINT2_DEFINED(eri,QD_x)
555  erieval->QD_x[v] = QDx;
556 #endif
557 #if LIBINT2_DEFINED(eri,QD_y)
558  erieval->QD_y[v] = QDy;
559 #endif
560 #if LIBINT2_DEFINED(eri,QD_z)
561  erieval->QD_z[v] = QDz;
562 #endif
563 
564 #if LIBINT2_DEFINED(eri,CD_x)
565  erieval->CD_x[v] = CD_x;
566 #endif
567 #if LIBINT2_DEFINED(eri,CD_y)
568  erieval->CD_y[v] = CD_y;
569 #endif
570 #if LIBINT2_DEFINED(eri,CD_z)
571  erieval->CD_z[v] = CD_z;
572 #endif
573 #if LIBINT2_DEFINED(eri,DC_x)
574  erieval->DC_x[v] = -CD_x;
575 #endif
576 #if LIBINT2_DEFINED(eri,DC_y)
577  erieval->DC_y[v] = -CD_y;
578 #endif
579 #if LIBINT2_DEFINED(eri,DC_z)
580  erieval->DC_z[v] = -CD_z;
581 #endif
582 #if LIBINT2_DEFINED(eri,oo2e)
583  erieval->oo2e[v] = 0.5*oogammaq;
584 #endif
585 
586  // Prefactors for interelectron transfer relation
587 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x)
588  erieval->TwoPRepITR_pfac0_0_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammap;
589 #endif
590 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y)
591  erieval->TwoPRepITR_pfac0_0_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammap;
592 #endif
593 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z)
594  erieval->TwoPRepITR_pfac0_0_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammap;
595 #endif
596 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x)
597  erieval->TwoPRepITR_pfac0_1_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammaq;
598 #endif
599 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y)
600  erieval->TwoPRepITR_pfac0_1_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammaq;
601 #endif
602 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z)
603  erieval->TwoPRepITR_pfac0_1_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammaq;
604 #endif
605 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x)
606  erieval->TwoPRepITR_pfac0_0_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammap;
607 #endif
608 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y)
609  erieval->TwoPRepITR_pfac0_0_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammap;
610 #endif
611 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z)
612  erieval->TwoPRepITR_pfac0_0_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammap;
613 #endif
614 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x)
615  erieval->TwoPRepITR_pfac0_1_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammaq;
616 #endif
617 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y)
618  erieval->TwoPRepITR_pfac0_1_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammaq;
619 #endif
620 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z)
621  erieval->TwoPRepITR_pfac0_1_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammaq;
622 #endif
623 #if LIBINT2_DEFINED(eri,eoz)
624  erieval->eoz[v] = gammaq*oogammap;
625 #endif
626 #if LIBINT2_DEFINED(eri,zoe)
627  erieval->zoe[v] = gammap*oogammaq;
628 #endif
629 
630  const double PQx = Px - Qx;
631  const double PQy = Py - Qy;
632  const double PQz = Pz - Qz;
633  const double PQ2 = PQx * PQx + PQy * PQy + PQz * PQz;
634  const double Wx = (gammap_o_gammapgammaq * Px + gammaq_o_gammapgammaq * Qx);
635  const double Wy = (gammap_o_gammapgammaq * Py + gammaq_o_gammapgammaq * Qy);
636  const double Wz = (gammap_o_gammapgammaq * Pz + gammaq_o_gammapgammaq * Qz);
637 
638 #if LIBINT2_DEFINED(eri,WP_x)
639  erieval->WP_x[v] = Wx - Px;
640 #endif
641 #if LIBINT2_DEFINED(eri,WP_y)
642  erieval->WP_y[v] = Wy - Py;
643 #endif
644 #if LIBINT2_DEFINED(eri,WP_z)
645  erieval->WP_z[v] = Wz - Pz;
646 #endif
647 #if LIBINT2_DEFINED(eri,WQ_x)
648  erieval->WQ_x[v] = Wx - Qx;
649 #endif
650 #if LIBINT2_DEFINED(eri,WQ_y)
651  erieval->WQ_y[v] = Wy - Qy;
652 #endif
653 #if LIBINT2_DEFINED(eri,WQ_z)
654  erieval->WQ_z[v] = Wz - Qz;
655 #endif
656 #if LIBINT2_DEFINED(eri,oo2ze)
657  erieval->oo2ze[v] = 0.5/(gammap+gammaq);
658 #endif
659 #if LIBINT2_DEFINED(eri,roz)
660  erieval->roz[v] = gammapq*oogammap;
661 #endif
662 #if LIBINT2_DEFINED(eri,roe)
663  erieval->roe[v] = gammapq*oogammaq;
664 #endif
665 
666  // prefactors for derivative ERI relations
667  if (deriv_order > 0) {
668 #if LIBINT2_DEFINED(eri,alpha1_rho_over_zeta2)
669  erieval->alpha1_rho_over_zeta2[v] = alpha0 * gammapq / (gammap * gammap);
670 #endif
671 #if LIBINT2_DEFINED(eri,alpha2_rho_over_zeta2)
672  erieval->alpha2_rho_over_zeta2[v] = alpha1 * gammapq / (gammap * gammap);
673 #endif
674 #if LIBINT2_DEFINED(eri,alpha3_rho_over_eta2)
675  erieval->alpha3_rho_over_eta2[v] = alpha2 * gammapq / (gammaq * gammaq);
676 #endif
677 #if LIBINT2_DEFINED(eri,alpha4_rho_over_eta2)
678  erieval->alpha4_rho_over_eta2[v] = alpha3 * gammapq / (gammaq * gammaq);
679 #endif
680 #if LIBINT2_DEFINED(eri,alpha1_over_zetapluseta)
681  erieval->alpha1_over_zetapluseta[v] = alpha0 / (gammap + gammaq);
682 #endif
683 #if LIBINT2_DEFINED(eri,alpha2_over_zetapluseta)
684  erieval->alpha2_over_zetapluseta[v] = alpha1 / (gammap + gammaq);
685 #endif
686 #if LIBINT2_DEFINED(eri,alpha3_over_zetapluseta)
687  erieval->alpha3_over_zetapluseta[v] = alpha2 / (gammap + gammaq);
688 #endif
689 #if LIBINT2_DEFINED(eri,alpha4_over_zetapluseta)
690  erieval->alpha4_over_zetapluseta[v] = alpha3 / (gammap + gammaq);
691 #endif
692 #if LIBINT2_DEFINED(eri,rho12_over_alpha1)
693  erieval->rho12_over_alpha1[v] = rhop / alpha0;
694 #endif
695 #if LIBINT2_DEFINED(eri,rho12_over_alpha2)
696  erieval->rho12_over_alpha2[v] = rhop / alpha1;
697 #endif
698 #if LIBINT2_DEFINED(eri,rho34_over_alpha3)
699  erieval->rho34_over_alpha3[v] = rhoq / alpha2;
700 #endif
701 #if LIBINT2_DEFINED(eri,rho34_over_alpha4)
702  erieval->rho34_over_alpha4[v] = rhoq / alpha3;
703 #endif
704 #if LIBINT2_DEFINED(eri,two_alpha0_bra)
705  erieval->two_alpha0_bra[v] = 2.0 * alpha0;
706 #endif
707 #if LIBINT2_DEFINED(eri,two_alpha0_ket)
708  erieval->two_alpha0_ket[v] = 2.0 * alpha1;
709 #endif
710 #if LIBINT2_DEFINED(eri,two_alpha1_bra)
711  erieval->two_alpha1_bra[v] = 2.0 * alpha2;
712 #endif
713 #if LIBINT2_DEFINED(eri,two_alpha1_ket)
714  erieval->two_alpha1_ket[v] = 2.0 * alpha3;
715 #endif
716  }
717 
718  const double K1 = exp(- rhop * AB2);
719  const double K2 = exp(- rhoq * CD2);
720  const double two_times_M_PI_to_25 = 34.986836655249725693;
721  double pfac = two_times_M_PI_to_25 * K1 * K2 / (gammap * gammaq * sqrt(gammap
722  + gammaq));
723  pfac *= c0 * c1 * c2 * c3;
724 
725  //calc_f(F, amtot, PQ2 * gammapq);
726  //libint2::FmEval_Reference2<double>::eval(F,PQ2*gammapq,amtot,1e-15);
727  fmeval_chebyshev.eval(F,PQ2*gammapq,amtot);
728  //fmeval_taylor.eval(F,PQ2*gammapq,amtot);
729 
730  // using dangerous macros from libint2.h
731 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(0))
732  erieval->LIBINT_T_SS_EREP_SS(0)[v] = pfac*F[0];
733 #endif
734 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(1))
735  erieval->LIBINT_T_SS_EREP_SS(1)[v] = pfac*F[1];
736 #endif
737 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(2))
738  erieval->LIBINT_T_SS_EREP_SS(2)[v] = pfac*F[2];
739 #endif
740 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(3))
741  erieval->LIBINT_T_SS_EREP_SS(3)[v] = pfac*F[3];
742 #endif
743 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(4))
744  erieval->LIBINT_T_SS_EREP_SS(4)[v] = pfac*F[4];
745 #endif
746 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(5))
747  erieval->LIBINT_T_SS_EREP_SS(5)[v] = pfac*F[5];
748 #endif
749 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(6))
750  erieval->LIBINT_T_SS_EREP_SS(6)[v] = pfac*F[6];
751 #endif
752 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(7))
753  erieval->LIBINT_T_SS_EREP_SS(7)[v] = pfac*F[7];
754 #endif
755 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(8))
756  erieval->LIBINT_T_SS_EREP_SS(8)[v] = pfac*F[8];
757 #endif
758 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(9))
759  erieval->LIBINT_T_SS_EREP_SS(9)[v] = pfac*F[9];
760 #endif
761 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(10))
762  erieval->LIBINT_T_SS_EREP_SS(10)[v] = pfac*F[10];
763 #endif
764 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(11))
765  erieval->LIBINT_T_SS_EREP_SS(11)[v] = pfac*F[11];
766 #endif
767 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(12))
768  erieval->LIBINT_T_SS_EREP_SS(12)[v] = pfac*F[12];
769 #endif
770 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(13))
771  erieval->LIBINT_T_SS_EREP_SS(13)[v] = pfac*F[13];
772 #endif
773 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(14))
774  erieval->LIBINT_T_SS_EREP_SS(14)[v] = pfac*F[14];
775 #endif
776 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(15))
777  erieval->LIBINT_T_SS_EREP_SS(15)[v] = pfac*F[15];
778 #endif
779 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(16))
780  erieval->LIBINT_T_SS_EREP_SS(16)[v] = pfac*F[16];
781 #endif
782 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(17))
783  erieval->LIBINT_T_SS_EREP_SS(17)[v] = pfac*F[17];
784 #endif
785 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(18))
786  erieval->LIBINT_T_SS_EREP_SS(18)[v] = pfac*F[18];
787 #endif
788 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(19))
789  erieval->LIBINT_T_SS_EREP_SS(19)[v] = pfac*F[19];
790 #endif
791 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(20))
792  erieval->LIBINT_T_SS_EREP_SS(20)[v] = pfac*F[20];
793 #endif
794  } // end of v loop
795 
796  }
797  }
798  } // end of primitive loops
799 }
800 
801 template<typename LibintEval>
802 void prep_libint2(LibintEval* erievals,
803  const RandomShellSet<2u>& rsqset,
804  int norm_flag,
805  int deriv_order = 0) {
806  const uint veclen = rsqset.exp[0].size();
807  const uint contrdepth = rsqset.exp[0][0].size();
808 
809  const unsigned int dummy_center1 = (LIBINT_SHELL_SET == LIBINT_SHELL_SET_STANDARD) ? 1 : 0;
810  const unsigned int dummy_center2 = (LIBINT_SHELL_SET == LIBINT_SHELL_SET_STANDARD) ? 3 : 2;
811  const double* A = &(rsqset.R[0][0]);
812  const double* B = &(rsqset.R[0][0]);
813  const double* C = &(rsqset.R[1][0]);
814  const double* D = &(rsqset.R[1][0]);
815 
816  const uint* am = rsqset.l;
817  const unsigned int amtot = am[0] + am[1] + deriv_order;
818  static double F[LIBINT_MAX_AM*4 + 6];
819 
820  uint p01 = 0;
821  for (uint p0 = 0; p0 < contrdepth; p0++) {
822  for (uint p1 = 0; p1 < contrdepth; p1++, p01++) {
823 
824  LibintEval* erieval = &erievals[p01];
825  erieval->veclen = veclen;
826 #if LIBINT2_FLOP_COUNT
827  erieval->nflops = erievals[0].nflops;
828 #endif
829 
830  for (uint v = 0; v < veclen; v++) {
831 
832  const double alpha0 = (dummy_center1 == 0) ? 0.0 : rsqset.exp[0][v][p0];
833  const double alpha1 = (dummy_center1 == 1) ? 0.0 : rsqset.exp[0][v][p0];
834  const double alpha2 = (dummy_center2 == 2) ? 0.0 : rsqset.exp[1][v][p1];
835  const double alpha3 = (dummy_center2 == 3) ? 0.0 : rsqset.exp[1][v][p1];
836 
837  const double c0 = (dummy_center1 == 0) ? 1.0 : rsqset.coef[0][v][p0];
838  const double c1 = (dummy_center1 == 1) ? 1.0 : rsqset.coef[0][v][p0];
839  const double c2 = (dummy_center2 == 2) ? 1.0 : rsqset.coef[1][v][p1];
840  const double c3 = (dummy_center2 == 3) ? 1.0 : rsqset.coef[1][v][p1];
841 
842  const double gammap = alpha0 + alpha1;
843  const double oogammap = 1.0 / gammap;
844  const double rhop = alpha0 * alpha1 * oogammap;
845  const double Px = (alpha0 * A[0] + alpha1 * B[0]) * oogammap;
846  const double Py = (alpha0 * A[1] + alpha1 * B[1]) * oogammap;
847  const double Pz = (alpha0 * A[2] + alpha1 * B[2]) * oogammap;
848  const double PAx = Px - A[0];
849  const double PAy = Py - A[1];
850  const double PAz = Pz - A[2];
851  const double PBx = Px - B[0];
852  const double PBy = Py - B[1];
853  const double PBz = Pz - B[2];
854  const double AB_x = A[0] - B[0];
855  const double AB_y = A[1] - B[1];
856  const double AB_z = A[2] - B[2];
857  const double AB2 = AB_x * AB_x + AB_y * AB_y + AB_z * AB_z;
858 
859 #if LIBINT2_DEFINED(eri,PA_x)
860  erieval->PA_x[v] = PAx;
861 #endif
862 #if LIBINT2_DEFINED(eri,PA_y)
863  erieval->PA_y[v] = PAy;
864 #endif
865 #if LIBINT2_DEFINED(eri,PA_z)
866  erieval->PA_z[v] = PAz;
867 #endif
868 #if LIBINT2_DEFINED(eri,PB_x)
869  erieval->PB_x[v] = PBx;
870 #endif
871 #if LIBINT2_DEFINED(eri,PB_y)
872  erieval->PB_y[v] = PBy;
873 #endif
874 #if LIBINT2_DEFINED(eri,PB_z)
875  erieval->PB_z[v] = PBz;
876 #endif
877 
878 #if LIBINT2_DEFINED(eri,AB_x)
879  erieval->AB_x[v] = AB_x;
880 #endif
881 #if LIBINT2_DEFINED(eri,AB_y)
882  erieval->AB_y[v] = AB_y;
883 #endif
884 #if LIBINT2_DEFINED(eri,AB_z)
885  erieval->AB_z[v] = AB_z;
886 #endif
887 #if LIBINT2_DEFINED(eri,oo2z)
888  erieval->oo2z[v] = 0.5*oogammap;
889 #endif
890 
891  const double gammaq = alpha2 + alpha3;
892  const double oogammaq = 1.0 / gammaq;
893  const double rhoq = alpha2 * alpha3 * oogammaq;
894  const double gammapq = gammap * gammaq / (gammap + gammaq);
895  const double gammap_o_gammapgammaq = gammapq * oogammaq;
896  const double gammaq_o_gammapgammaq = gammapq * oogammap;
897  const double Qx = (alpha2 * C[0] + alpha3 * D[0]) * oogammaq;
898  const double Qy = (alpha2 * C[1] + alpha3 * D[1]) * oogammaq;
899  const double Qz = (alpha2 * C[2] + alpha3 * D[2]) * oogammaq;
900  const double QCx = Qx - C[0];
901  const double QCy = Qy - C[1];
902  const double QCz = Qz - C[2];
903  const double QDx = Qx - D[0];
904  const double QDy = Qy - D[1];
905  const double QDz = Qz - D[2];
906  const double CD_x = C[0] - D[0];
907  const double CD_y = C[1] - D[1];
908  const double CD_z = C[2] - D[2];
909  const double CD2 = CD_x * CD_x + CD_y * CD_y + CD_z * CD_z;
910 
911 #if LIBINT2_DEFINED(eri,QC_x)
912  erieval->QC_x[v] = QCx;
913 #endif
914 #if LIBINT2_DEFINED(eri,QC_y)
915  erieval->QC_y[v] = QCy;
916 #endif
917 #if LIBINT2_DEFINED(eri,QC_z)
918  erieval->QC_z[v] = QCz;
919 #endif
920 #if LIBINT2_DEFINED(eri,QD_x)
921  erieval->QD_x[v] = QDx;
922 #endif
923 #if LIBINT2_DEFINED(eri,QD_y)
924  erieval->QD_y[v] = QDy;
925 #endif
926 #if LIBINT2_DEFINED(eri,QD_z)
927  erieval->QD_z[v] = QDz;
928 #endif
929 
930 #if LIBINT2_DEFINED(eri,CD_x)
931  erieval->CD_x[v] = CD_x;
932 #endif
933 #if LIBINT2_DEFINED(eri,CD_y)
934  erieval->CD_y[v] = CD_y;
935 #endif
936 #if LIBINT2_DEFINED(eri,CD_z)
937  erieval->CD_z[v] = CD_z;
938 #endif
939 #if LIBINT2_DEFINED(eri,oo2e)
940  erieval->oo2e[v] = 0.5*oogammaq;
941 #endif
942 
943  // Prefactors for interelectron transfer relation
944 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x)
945  erieval->TwoPRepITR_pfac0_0_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammap;
946 #endif
947 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y)
948  erieval->TwoPRepITR_pfac0_0_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammap;
949 #endif
950 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z)
951  erieval->TwoPRepITR_pfac0_0_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammap;
952 #endif
953 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x)
954  erieval->TwoPRepITR_pfac0_1_0_x[v] = - (alpha1*AB_x + alpha3*CD_x)*oogammaq;
955 #endif
956 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y)
957  erieval->TwoPRepITR_pfac0_1_0_y[v] = - (alpha1*AB_y + alpha3*CD_y)*oogammaq;
958 #endif
959 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z)
960  erieval->TwoPRepITR_pfac0_1_0_z[v] = - (alpha1*AB_z + alpha3*CD_z)*oogammaq;
961 #endif
962 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x)
963  erieval->TwoPRepITR_pfac0_0_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammap;
964 #endif
965 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y)
966  erieval->TwoPRepITR_pfac0_0_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammap;
967 #endif
968 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z)
969  erieval->TwoPRepITR_pfac0_0_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammap;
970 #endif
971 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x)
972  erieval->TwoPRepITR_pfac0_1_1_x[v] = (alpha0*AB_x + alpha2*CD_x)*oogammaq;
973 #endif
974 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y)
975  erieval->TwoPRepITR_pfac0_1_1_y[v] = (alpha0*AB_y + alpha2*CD_y)*oogammaq;
976 #endif
977 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z)
978  erieval->TwoPRepITR_pfac0_1_1_z[v] = (alpha0*AB_z + alpha2*CD_z)*oogammaq;
979 #endif
980 #if LIBINT2_DEFINED(eri,eoz)
981  erieval->eoz[v] = gammaq*oogammap;
982 #endif
983 #if LIBINT2_DEFINED(eri,zoe)
984  erieval->zoe[v] = gammap*oogammaq;
985 #endif
986 
987  const double PQx = Px - Qx;
988  const double PQy = Py - Qy;
989  const double PQz = Pz - Qz;
990  const double PQ2 = PQx * PQx + PQy * PQy + PQz * PQz;
991  const double Wx = (gammap_o_gammapgammaq * Px + gammaq_o_gammapgammaq * Qx);
992  const double Wy = (gammap_o_gammapgammaq * Py + gammaq_o_gammapgammaq * Qy);
993  const double Wz = (gammap_o_gammapgammaq * Pz + gammaq_o_gammapgammaq * Qz);
994 
995 #if LIBINT2_DEFINED(eri,WP_x)
996  erieval->WP_x[v] = Wx - Px;
997 #endif
998 #if LIBINT2_DEFINED(eri,WP_y)
999  erieval->WP_y[v] = Wy - Py;
1000 #endif
1001 #if LIBINT2_DEFINED(eri,WP_z)
1002  erieval->WP_z[v] = Wz - Pz;
1003 #endif
1004 #if LIBINT2_DEFINED(eri,WQ_x)
1005  erieval->WQ_x[v] = Wx - Qx;
1006 #endif
1007 #if LIBINT2_DEFINED(eri,WQ_y)
1008  erieval->WQ_y[v] = Wy - Qy;
1009 #endif
1010 #if LIBINT2_DEFINED(eri,WQ_z)
1011  erieval->WQ_z[v] = Wz - Qz;
1012 #endif
1013 #if LIBINT2_DEFINED(eri,oo2ze)
1014  erieval->oo2ze[v] = 0.5/(gammap+gammaq);
1015 #endif
1016 #if LIBINT2_DEFINED(eri,roz)
1017  erieval->roz[v] = gammapq*oogammap;
1018 #endif
1019 #if LIBINT2_DEFINED(eri,roe)
1020  erieval->roe[v] = gammapq*oogammaq;
1021 #endif
1022 
1023  // prefactors for derivative ERI relations
1024  if (deriv_order > 0) {
1025 #if LIBINT2_DEFINED(eri,alpha1_rho_over_zeta2)
1026  erieval->alpha1_rho_over_zeta2[v] = alpha0 * gammapq / (gammap * gammap);
1027 #endif
1028 #if LIBINT2_DEFINED(eri,alpha2_rho_over_zeta2)
1029  erieval->alpha2_rho_over_zeta2[v] = alpha1 * gammapq / (gammap * gammap);
1030 #endif
1031 #if LIBINT2_DEFINED(eri,alpha3_rho_over_eta2)
1032  erieval->alpha3_rho_over_eta2[v] = alpha2 * gammapq / (gammaq * gammaq);
1033 #endif
1034 #if LIBINT2_DEFINED(eri,alpha4_rho_over_eta2)
1035  erieval->alpha4_rho_over_eta2[v] = alpha3 * gammapq / (gammaq * gammaq);
1036 #endif
1037 #if LIBINT2_DEFINED(eri,alpha1_over_zetapluseta)
1038  erieval->alpha1_over_zetapluseta[v] = alpha0 / (gammap + gammaq);
1039 #endif
1040 #if LIBINT2_DEFINED(eri,alpha2_over_zetapluseta)
1041  erieval->alpha2_over_zetapluseta[v] = alpha1 / (gammap + gammaq);
1042 #endif
1043 #if LIBINT2_DEFINED(eri,alpha3_over_zetapluseta)
1044  erieval->alpha3_over_zetapluseta[v] = alpha2 / (gammap + gammaq);
1045 #endif
1046 #if LIBINT2_DEFINED(eri,alpha4_over_zetapluseta)
1047  erieval->alpha4_over_zetapluseta[v] = alpha3 / (gammap + gammaq);
1048 #endif
1049 #if LIBINT2_DEFINED(eri,rho12_over_alpha1)
1050  erieval->rho12_over_alpha1[v] = rhop / alpha0;
1051 #endif
1052 #if LIBINT2_DEFINED(eri,rho12_over_alpha2)
1053  erieval->rho12_over_alpha2[v] = rhop / alpha1;
1054 #endif
1055 #if LIBINT2_DEFINED(eri,rho34_over_alpha3)
1056  erieval->rho34_over_alpha3[v] = rhoq / alpha2;
1057 #endif
1058 #if LIBINT2_DEFINED(eri,rho34_over_alpha4)
1059  erieval->rho34_over_alpha4[v] = rhoq / alpha3;
1060 #endif
1061 #if LIBINT2_DEFINED(eri,two_alpha0_bra)
1062  erieval->two_alpha0_bra[v] = 2.0 * alpha0;
1063 #endif
1064 #if LIBINT2_DEFINED(eri,two_alpha0_ket)
1065  erieval->two_alpha0_ket[v] = 2.0 * alpha1;
1066 #endif
1067 #if LIBINT2_DEFINED(eri,two_alpha1_bra)
1068  erieval->two_alpha1_bra[v] = 2.0 * alpha2;
1069 #endif
1070 #if LIBINT2_DEFINED(eri,two_alpha1_ket)
1071  erieval->two_alpha1_ket[v] = 2.0 * alpha3;
1072 #endif
1073  }
1074 
1075  const double K1 = exp(- rhop * AB2 );
1076  const double K2 = exp(- rhoq * CD2 );
1077  const double two_times_M_PI_to_25 = 34.986836655249725693;
1078  double pfac = two_times_M_PI_to_25 * K1 * K2 / (gammap * gammaq * sqrt(gammap
1079  + gammaq));
1080  pfac *= c0 * c1 * c2 * c3;
1081 
1082  //calc_f(F, amtot, PQ2 * gammapq);
1083  //libint2::FmEval_Reference2<double>::eval(F,PQ2*gammapq,amtot,1e-15);
1084  fmeval_chebyshev.eval(F,PQ2*gammapq,amtot);
1085  //fmeval_taylor.eval(F,PQ2*gammapq,amtot);
1086 
1087  // using dangerous macros from libint2.h
1088 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(0))
1089  erieval->LIBINT_T_SS_EREP_SS(0)[v] = pfac*F[0];
1090 #endif
1091 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(1))
1092  erieval->LIBINT_T_SS_EREP_SS(1)[v] = pfac*F[1];
1093 #endif
1094 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(2))
1095  erieval->LIBINT_T_SS_EREP_SS(2)[v] = pfac*F[2];
1096 #endif
1097 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(3))
1098  erieval->LIBINT_T_SS_EREP_SS(3)[v] = pfac*F[3];
1099 #endif
1100 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(4))
1101  erieval->LIBINT_T_SS_EREP_SS(4)[v] = pfac*F[4];
1102 #endif
1103 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(5))
1104  erieval->LIBINT_T_SS_EREP_SS(5)[v] = pfac*F[5];
1105 #endif
1106 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(6))
1107  erieval->LIBINT_T_SS_EREP_SS(6)[v] = pfac*F[6];
1108 #endif
1109 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(7))
1110  erieval->LIBINT_T_SS_EREP_SS(7)[v] = pfac*F[7];
1111 #endif
1112 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(8))
1113  erieval->LIBINT_T_SS_EREP_SS(8)[v] = pfac*F[8];
1114 #endif
1115 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(9))
1116  erieval->LIBINT_T_SS_EREP_SS(9)[v] = pfac*F[9];
1117 #endif
1118 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(10))
1119  erieval->LIBINT_T_SS_EREP_SS(10)[v] = pfac*F[10];
1120 #endif
1121 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(11))
1122  erieval->LIBINT_T_SS_EREP_SS(11)[v] = pfac*F[11];
1123 #endif
1124 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(12))
1125  erieval->LIBINT_T_SS_EREP_SS(12)[v] = pfac*F[12];
1126 #endif
1127 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(13))
1128  erieval->LIBINT_T_SS_EREP_SS(13)[v] = pfac*F[13];
1129 #endif
1130 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(14))
1131  erieval->LIBINT_T_SS_EREP_SS(14)[v] = pfac*F[14];
1132 #endif
1133 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(15))
1134  erieval->LIBINT_T_SS_EREP_SS(15)[v] = pfac*F[15];
1135 #endif
1136 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(16))
1137  erieval->LIBINT_T_SS_EREP_SS(16)[v] = pfac*F[16];
1138 #endif
1139 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(17))
1140  erieval->LIBINT_T_SS_EREP_SS(17)[v] = pfac*F[17];
1141 #endif
1142 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(18))
1143  erieval->LIBINT_T_SS_EREP_SS(18)[v] = pfac*F[18];
1144 #endif
1145 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(19))
1146  erieval->LIBINT_T_SS_EREP_SS(19)[v] = pfac*F[19];
1147 #endif
1148 #if LIBINT2_DEFINED(eri,LIBINT_T_SS_EREP_SS(20))
1149  erieval->LIBINT_T_SS_EREP_SS(20)[v] = pfac*F[20];
1150 #endif
1151  } // end of v loop
1152 
1153  }
1154  } // end of primitive loops
1155 }
1156 
Definition: prep_libint2.h:54
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
void eval(Real *Fm, Real x, int m_max) const
fills in Fm with computed Boys function values for m in [0,mmax]
Definition: boys.h:295
Computes the Boys function, $ F_m (T) = ^1 u^{2m} (-T u^2) \, { d}u $, using Taylor interpolation of ...
Definition: boys.h:752
Definition: prep_libint2.h:35