20 #ifndef _libint2_src_bin_testeri_eri_h_ 21 #define _libint2_src_bin_testeri_eri_h_ 30 #include <libint2/boys.h> 32 #if defined(__cplusplus) 37 typedef mpf_class LIBINT2_REF_REALTYPE;
39 inline mpf_class exp(mpf_class x) {
40 mpfr_t x_r; mpfr_init(x_r);
41 mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
43 mpfr_t expx_r; mpfr_init(expx_r);
44 mpfr_exp(expx_r, x_r, MPFR_RNDN);
48 mpfr_get_f(expx, expx_r, MPFR_RNDN);
49 mpf_class result(expx);
53 inline mpf_class pow(mpf_class x,
int a) {
57 mpf_pow_ui(x_to_a, x.get_mpf_t(), (
unsigned int) a);
59 mpf_pow_ui(x_to_a, x.get_mpf_t(), (
unsigned int) (-a));
60 mpf_class result(x_to_a);
62 result = 1.0 / result;
66 inline mpf_class erf(mpf_class x) {
67 mpfr_t x_r; mpfr_init(x_r);
68 mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
70 mpfr_t erfx_r; mpfr_init(erfx_r);
71 mpfr_erf(erfx_r, x_r, MPFR_RNDN);
75 mpfr_get_f(erfx, erfx_r, MPFR_RNDN);
76 mpf_class result(erfx);
80 inline mpf_class log(mpf_class x) {
81 mpfr_t x_r; mpfr_init(x_r);
82 mpfr_set_f(x_r, x.get_mpf_t(), MPFR_RNDN);
84 mpfr_t logx_r; mpfr_init(logx_r);
85 mpfr_log(logx_r, x_r, MPFR_RNDN);
89 mpfr_get_f(logx, logx_r, MPFR_RNDN);
90 mpf_class result(logx);
95 typedef double LIBINT2_REF_REALTYPE;
98 typedef double LIBINT2_REF_REALTYPE;
105 #if defined(__clusplus) 106 typedef mpz_class BigInt;
108 typedef double BigInt;
113 const static unsigned int MAXFAC = 100;
116 df =
new BigInt[2*MAXFAC];
120 for(
unsigned long i=3; i<MAXFAC*2; i++) {
121 df[i] = (i-1)*df[i-2];
124 fac =
new BigInt[MAXFAC];
126 for(
unsigned long i=1;i<MAXFAC;i++)
129 bc =
new BigInt*[MAXFAC];
130 bc[0] =
new BigInt[MAXFAC*MAXFAC];
131 for(
unsigned long i=1;i<MAXFAC;i++)
132 bc[i] = bc[i-1] + MAXFAC;
133 for(
unsigned long i=0;i<MAXFAC;i++) {
135 for(
unsigned long j=1;j<=i;j++) {
136 bc[i][j] = (bc[i][j-1] * (i-j+1) )/ j;
148 template <
typename Real>
149 Real norm_const(
unsigned int l1,
unsigned int m1,
unsigned int n1,
150 Real alpha1,
const Real* A)
152 const Real sqrt_twoalpha1_over_pi = sqrt(2*alpha1/M_PI);
153 const Real sqrtsqrt_twoalpha1_over_pi = sqrt(sqrt_twoalpha1_over_pi);
154 const Real sqrt_alpha1 = sqrt(alpha1);
156 return sqrt_twoalpha1_over_pi * sqrtsqrt_twoalpha1_over_pi * (2 * pow(sqrt_alpha1,l1+m1+n1))/sqrt(df[2*l1]*df[2*m1]*df[2*n1]);
170 template <
typename Real>
void calc_f(Real* F,
int n, Real t)
179 static Real K = 1.0/M_2_SQRTPI;
187 for(m=0; m<=n-1; m++){
188 F[m+1] = ((2*m + 1)*F[m] - et)/(t2);
198 num = expmath.df[m2];
204 term1 = num/expmath.df[m2+2*i+2];
206 }
while (term1 > EPS && i < expmath.MAXFAC);
209 F[m] = (t2*F[m+1] + et)/(2*m+1);
215 template <
typename Int>
217 if (a%2 == 1)
return -1;
223 int min(
int t1,
unsigned int t2) {
return t1 > (int)t2 ? (
int)t2 : t1; }
236 inline LIBINT2_REF_REALTYPE eri(
unsigned int l1,
unsigned int m1,
unsigned int n1,
237 LIBINT2_REF_REALTYPE alpha1,
const LIBINT2_REF_REALTYPE* A,
238 unsigned int l2,
unsigned int m2,
unsigned int n2,
239 LIBINT2_REF_REALTYPE alpha2,
const LIBINT2_REF_REALTYPE* B,
240 unsigned int l3,
unsigned int m3,
unsigned int n3,
241 LIBINT2_REF_REALTYPE alpha3,
const LIBINT2_REF_REALTYPE* C,
242 unsigned int l4,
unsigned int m4,
unsigned int n4,
243 LIBINT2_REF_REALTYPE alpha4,
const LIBINT2_REF_REALTYPE* D,
248 const LIBINT2_REF_REALTYPE gammap = alpha1 + alpha2;
249 const LIBINT2_REF_REALTYPE Px = (alpha1*A[0] + alpha2*B[0])/gammap;
250 const LIBINT2_REF_REALTYPE Py = (alpha1*A[1] + alpha2*B[1])/gammap;
251 const LIBINT2_REF_REALTYPE Pz = (alpha1*A[2] + alpha2*B[2])/gammap;
252 const LIBINT2_REF_REALTYPE PAx = Px - A[0];
253 const LIBINT2_REF_REALTYPE PAy = Py - A[1];
254 const LIBINT2_REF_REALTYPE PAz = Pz - A[2];
255 const LIBINT2_REF_REALTYPE PBx = Px - B[0];
256 const LIBINT2_REF_REALTYPE PBy = Py - B[1];
257 const LIBINT2_REF_REALTYPE PBz = Pz - B[2];
258 const LIBINT2_REF_REALTYPE AB2 = (A[0]-B[0])*(A[0]-B[0]) + (A[1]-B[1])*(A[1]-B[1]) + (A[2]-B[2])*(A[2]-B[2]);
260 const LIBINT2_REF_REALTYPE gammaq = alpha3 + alpha4;
261 const LIBINT2_REF_REALTYPE gammapq = gammap*gammaq/(gammap+gammaq);
262 const LIBINT2_REF_REALTYPE Qx = (alpha3*C[0] + alpha4*D[0])/gammaq;
263 const LIBINT2_REF_REALTYPE Qy = (alpha3*C[1] + alpha4*D[1])/gammaq;
264 const LIBINT2_REF_REALTYPE Qz = (alpha3*C[2] + alpha4*D[2])/gammaq;
265 const LIBINT2_REF_REALTYPE QCx = Qx - C[0];
266 const LIBINT2_REF_REALTYPE QCy = Qy - C[1];
267 const LIBINT2_REF_REALTYPE QCz = Qz - C[2];
268 const LIBINT2_REF_REALTYPE QDx = Qx - D[0];
269 const LIBINT2_REF_REALTYPE QDy = Qy - D[1];
270 const LIBINT2_REF_REALTYPE QDz = Qz - D[2];
271 const LIBINT2_REF_REALTYPE CD2 = (C[0]-D[0])*(C[0]-D[0]) + (C[1]-D[1])*(C[1]-D[1]) + (C[2]-D[2])*(C[2]-D[2]);
273 const LIBINT2_REF_REALTYPE PQx = Px - Qx;
274 const LIBINT2_REF_REALTYPE PQy = Py - Qy;
275 const LIBINT2_REF_REALTYPE PQz = Pz - Qz;
276 const LIBINT2_REF_REALTYPE PQ2 = PQx*PQx + PQy*PQy + PQz*PQz;
278 const LIBINT2_REF_REALTYPE four = 4;
280 int u1,u2,v1,v2,w1,w2,tx,ty,tz,txmax,tymax,tzmax;
282 int lp,lq,mp,mq,np,nq;
284 LIBINT2_REF_REALTYPE *flp, *flq, *fmp, *fmq, *fnp, *fnq;
285 LIBINT2_REF_REALTYPE *F;
286 LIBINT2_REF_REALTYPE K1, K2;
287 LIBINT2_REF_REALTYPE Gx,Gy,Gz;
288 LIBINT2_REF_REALTYPE pfac;
289 LIBINT2_REF_REALTYPE result = 0.0;
290 LIBINT2_REF_REALTYPE tmp;
291 int u1max,u2max,v1max,v2max,w1max,w2max;
293 K1 = exp(-alpha1*alpha2*AB2/gammap);
294 K2 = exp(-alpha3*alpha4*CD2/gammaq);
295 const LIBINT2_REF_REALTYPE m_pi = M_PI;
296 const LIBINT2_REF_REALTYPE m_pi_2 = m_pi * m_pi;
297 pfac = 2*sqrt(m_pi_2 * m_pi_2 * m_pi)*K1*K2/(gammap*gammaq*sqrt(gammap+gammaq));
300 pfac *= expmath.norm_const(l1,m1,n1,alpha1,A);
301 pfac *= expmath.norm_const(l2,m2,n2,alpha2,B);
302 pfac *= expmath.norm_const(l3,m3,n3,alpha3,C);
303 pfac *= expmath.norm_const(l4,m4,n4,alpha4,D);
306 const int ltot = l1+l2+l3+l4+m1+m2+m3+m4+n1+n2+n3+n4;
307 F =
new LIBINT2_REF_REALTYPE[ltot+1];
313 std::cout <<
"pfac = " << pfac <<
" K1 = " << K1 <<
" K2 = " << K2 << std::endl;
314 std::cout <<
"PQ2 = " << PQ2 <<
" gammapq = " << gammapq << std::endl;
315 for(
int i=0; i<=ltot; ++i)
316 std::cout <<
"F[" << i <<
"] = " << F[i] << std::endl;
319 flp =
new LIBINT2_REF_REALTYPE[l1+l2+1];
320 for(k=0;k<=l1+l2;k++) {
322 for(i=0;i<=libint2::min(k,l1);i++) {
324 if (j > l2)
continue;
325 tmp = expmath.bc[l1][i]*expmath.bc[l2][j];
327 tmp *= pow(PAx, l1 - i);
329 tmp *= pow(PBx, l2 - j);
333 fmp =
new LIBINT2_REF_REALTYPE[m1+m2+1];
334 for(k=0;k<=m1+m2;k++) {
336 for(i=0;i<=libint2::min(k,m1);i++) {
338 if (j > m2)
continue;
339 tmp = expmath.bc[m1][i]*expmath.bc[m2][j];
341 tmp *= pow(PAy, m1 - i);
343 tmp *= pow(PBy, m2 - j);
347 fnp =
new LIBINT2_REF_REALTYPE[n1+n2+1];
348 for(k=0;k<=n1+n2;k++) {
350 for(i=0;i<=libint2::min(k,n1);i++) {
352 if (j > n2)
continue;
353 tmp = expmath.bc[n1][i]*expmath.bc[n2][j];
355 tmp *= pow(PAz, n1 - i);
357 tmp *= pow(PBz, n2 - j);
361 flq =
new LIBINT2_REF_REALTYPE[l3+l4+1];
362 for(k=0;k<=l3+l4;k++) {
364 for(i=0;i<=libint2::min(k,l3);i++) {
366 if (j > l4)
continue;
367 tmp = expmath.bc[l3][i]*expmath.bc[l4][j];
369 tmp *= pow(QCx, l3 - i);
371 tmp *= pow(QDx, l4 - j);
375 fmq =
new LIBINT2_REF_REALTYPE[m3+m4+1];
376 for(k=0;k<=m3+m4;k++) {
378 for(i=0;i<=libint2::min(k,m3);i++) {
380 if (j > m4)
continue;
381 tmp = expmath.bc[m3][i]*expmath.bc[m4][j];
383 tmp *= pow(QCy, m3 - i);
385 tmp *= pow(QDy, m4 - j);
389 fnq =
new LIBINT2_REF_REALTYPE[n3+n4+1];
390 for(k=0;k<=n3+n4;k++) {
392 for(i=0;i<=libint2:: min(k,n3);i++) {
394 if (j > n4)
continue;
395 tmp = expmath.bc[n3][i]*expmath.bc[n4][j];
397 tmp *= pow(QCz, n3 - i);
399 tmp *= pow(QDz, n4 - j);
404 for (lp = 0; lp <= l1 + l2; lp++) {
405 for (lq = 0; lq <= l3 + l4; lq++) {
408 for (u1 = 0; u1 <= u1max; u1++) {
409 for (u2 = 0; u2 <= u2max; u2++) {
410 Gx = parity(lp) * flp[lp] * flq[lq] * expmath.fac[lp]
411 * expmath.fac[lq] * pow(gammap, u1 - lp) * pow(gammaq, u2 - lq)
412 * expmath.fac[lp + lq - 2 * u1 - 2 * u2]
413 * pow(gammapq, lp + lq - 2 * u1 - 2 * u2)
414 / (expmath.fac[u1] * expmath.fac[u2] * expmath.fac[lp - 2 * u1]
415 * expmath.fac[lq - 2 * u2]);
416 for (mp = 0; mp <= m1 + m2; mp++) {
417 for (mq = 0; mq <= m3 + m4; mq++) {
420 for (v1 = 0; v1 <= v1max; v1++) {
421 for (v2 = 0; v2 <= v2max; v2++) {
423 parity(mp) * fmp[mp] * fmq[mq] * expmath.fac[mp]
424 * expmath.fac[mq] * pow(gammap, v1 - mp)
425 * pow(gammaq, v2 - mq)
426 * expmath.fac[mp + mq - 2 * v1 - 2 * v2]
427 * pow(gammapq, mp + mq - 2 * v1 - 2 * v2)
428 / (expmath.fac[v1] * expmath.fac[v2]
429 * expmath.fac[mp - 2 * v1]
430 * expmath.fac[mq - 2 * v2]);
431 for (np = 0; np <= n1 + n2; np++) {
432 for (nq = 0; nq <= n3 + n4; nq++) {
435 for (w1 = 0; w1 <= w1max; w1++) {
436 for (w2 = 0; w2 <= w2max; w2++) {
437 Gz = parity(np) * fnp[np] * fnq[nq] * expmath.fac[np]
438 * expmath.fac[nq] * pow(gammap, w1 - np)
439 * pow(gammaq, w2 - nq)
440 * expmath.fac[np + nq - 2 * w1 - 2 * w2]
441 * pow(gammapq, np + nq - 2 * w1 - 2 * w2)
442 / (expmath.fac[w1] * expmath.fac[w2]
443 * expmath.fac[np - 2 * w1]
444 * expmath.fac[nq - 2 * w2]);
445 txmax = (lp + lq - 2 * u1 - 2 * u2) / 2;
446 tymax = (mp + mq - 2 * v1 - 2 * v2) / 2;
447 tzmax = (np + nq - 2 * w1 - 2 * w2) / 2;
448 for (tx = 0; tx <= txmax; tx++) {
449 for (ty = 0; ty <= tymax; ty++) {
450 for (tz = 0; tz <= tzmax; tz++) {
451 zeta = lp + lq + mp + mq + np + nq - 2 * u1
452 - 2 * u2 - 2 * v1 - 2 * v2 - 2 * w1 - 2 * w2
454 result += Gx * Gy * Gz * F[zeta]
455 * parity(tx + ty + tz)
457 lp + lq - 2 * u1 - 2 * u2 - 2 * tx)
459 mp + mq - 2 * v1 - 2 * v2 - 2 * ty)
461 np + nq - 2 * w1 - 2 * w2 - 2 * tz)
464 u1 + u2 + tx + v1 + v2 + ty + w1 + w2
465 + tz) * pow(gammapq, tx)
466 * pow(gammapq, ty) * pow(gammapq, tz)
467 * expmath.fac[lp + lq - 2 * u1 - 2 * u2
468 - 2 * tx] * expmath.fac[tx]
469 * expmath.fac[mp + mq - 2 * v1 - 2 * v2
470 - 2 * ty] * expmath.fac[ty]
471 * expmath.fac[np + nq - 2 * w1 - 2 * w2
472 - 2 * tz] * expmath.fac[tz]);
474 std::cout <<
"Gx = " << Gx << std::endl;
475 std::cout <<
"Gy = " << Gy << std::endl;
476 std::cout <<
"Gz = " << Gz << std::endl;
477 std::cout <<
"result = " << result << std::endl;
510 inline LIBINT2_REF_REALTYPE eri(
const unsigned int* deriv_index,
511 unsigned int l1,
unsigned int m1,
unsigned int n1,
512 LIBINT2_REF_REALTYPE alpha1,
const LIBINT2_REF_REALTYPE* A,
513 unsigned int l2,
unsigned int m2,
unsigned int n2,
514 LIBINT2_REF_REALTYPE alpha2,
const LIBINT2_REF_REALTYPE* B,
515 unsigned int l3,
unsigned int m3,
unsigned int n3,
516 LIBINT2_REF_REALTYPE alpha3,
const LIBINT2_REF_REALTYPE* C,
517 unsigned int l4,
unsigned int m4,
unsigned int n4,
518 LIBINT2_REF_REALTYPE alpha4,
const LIBINT2_REF_REALTYPE* D,
520 const unsigned int deriv_order = std::accumulate(deriv_index, deriv_index+12, 0u);
521 if (deriv_order == 0)
522 return eri(l1, m1, n1, alpha1, A,
523 l2, m2, n2, alpha2, B,
524 l3, m3, n3, alpha3, C,
525 l4, m4, n4, alpha4, D,
530 bool operator()(
unsigned int val) {
537 const unsigned int* di_ptr = std::find_if(deriv_index, deriv_index+12, nonzero);
538 const unsigned int di = di_ptr - deriv_index;
539 const unsigned int dc = di / 3;
540 const unsigned int dxyz = di % 3;
542 unsigned int qn[4][3];
556 LIBINT2_REF_REALTYPE alpha[4];
563 unsigned int deriv_index_minus1[12]; std::copy(deriv_index, deriv_index+12, deriv_index_minus1);
564 --deriv_index_minus1[di];
565 LIBINT2_REF_REALTYPE result = eri(deriv_index_minus1,
566 qn[0][0], qn[0][1], qn[0][2], alpha[0], A,
567 qn[1][0], qn[1][1], qn[1][2], alpha[1], B,
568 qn[2][0], qn[2][1], qn[2][2], alpha[2], C,
569 qn[3][0], qn[3][1], qn[3][2], alpha[3], D,
574 if (qn[dc][dxyz] > 0) {
576 result -= eri(deriv_index_minus1,
577 qn[0][0], qn[0][1], qn[0][2], alpha[0], A,
578 qn[1][0], qn[1][1], qn[1][2], alpha[1], B,
579 qn[2][0], qn[2][1], qn[2][2], alpha[2], C,
580 qn[3][0], qn[3][1], qn[3][2], alpha[3], D,
593 #endif // header guard
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
static void eval(Real *Fm, Real t, size_t mmax, Real absolute_precision)
fills up an array of Fm(T) for m in [0,mmax]
Definition: boys.h:212