LIBINT  2.1.0-stable
eri.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 _libint2_src_bin_testeri_eri_h_
21 #define _libint2_src_bin_testeri_eri_h_
22 
23 #include <cmath>
24 #include <numeric>
25 #include <algorithm>
26 #include <cassert>
27 #include <iostream>
28 #include <stdlib.h>
29 
30 #include <libint2/boys.h>
31 
32 #if defined(__cplusplus)
33 #if HAVE_MPFR
34 # include <cstddef>
35 # include <gmpxx.h>
36 # include <mpfr.h>
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);
42 
43  mpfr_t expx_r; mpfr_init(expx_r);
44  mpfr_exp(expx_r, x_r, MPFR_RNDN);
45 
46  mpf_t expx;
47  mpf_init(expx);
48  mpfr_get_f(expx, expx_r, MPFR_RNDN);
49  mpf_class result(expx);
50  return result;
51  }
53  inline mpf_class pow(mpf_class x, int a) {
54  mpf_t x_to_a;
55  mpf_init(x_to_a);
56  if (a >= 0)
57  mpf_pow_ui(x_to_a, x.get_mpf_t(), (unsigned int) a);
58  else
59  mpf_pow_ui(x_to_a, x.get_mpf_t(), (unsigned int) (-a));
60  mpf_class result(x_to_a);
61  if (a < 0)
62  result = 1.0 / result;
63  return result;
64  }
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);
69 
70  mpfr_t erfx_r; mpfr_init(erfx_r);
71  mpfr_erf(erfx_r, x_r, MPFR_RNDN);
72 
73  mpf_t erfx;
74  mpf_init(erfx);
75  mpfr_get_f(erfx, erfx_r, MPFR_RNDN);
76  mpf_class result(erfx);
77  return result;
78  }
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);
83 
84  mpfr_t logx_r; mpfr_init(logx_r);
85  mpfr_log(logx_r, x_r, MPFR_RNDN);
86 
87  mpf_t logx;
88  mpf_init(logx);
89  mpfr_get_f(logx, logx_r, MPFR_RNDN);
90  mpf_class result(logx);
91  return result;
92  }
93 
94 #else
95  typedef double LIBINT2_REF_REALTYPE;
96 #endif
97 #else
98  typedef double LIBINT2_REF_REALTYPE;
99 #endif
100 
101 #define EPS 1.0E-17 /* Absolute precision in computing Fm(t)
102  (see recursion:calc_fij() ) */
103 
105 #if defined(__clusplus)
106  typedef mpz_class BigInt;
107 #else
108  typedef double BigInt;
109 #endif
110  BigInt *df;
111  BigInt *fac;
112  BigInt **bc;
113  const static unsigned int MAXFAC = 100;
114 
115  ExpensiveMath() {
116  df = new BigInt[2*MAXFAC];
117  df[0] = 1ul;
118  df[1] = 1ul;
119  df[2] = 1ul;
120  for(unsigned long i=3; i<MAXFAC*2; i++) {
121  df[i] = (i-1)*df[i-2];
122  }
123 
124  fac = new BigInt[MAXFAC];
125  fac[0] = 1ul;
126  for(unsigned long i=1;i<MAXFAC;i++)
127  fac[i] = fac[i-1]*i;
128 
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++) {
134  bc[i][0] = 1ul;
135  for(unsigned long j=1;j<=i;j++) {
136  bc[i][j] = (bc[i][j-1] * (i-j+1) )/ j; // fac[i]/(fac[i-j]*fac[j]);
137  }
138  }
139  }
140 
141  ~ExpensiveMath() {
142  delete[] df;
143  delete[] fac;
144  delete[] bc[0];
145  delete[] bc;
146  }
147 
148  template <typename Real>
149  Real norm_const(unsigned int l1, unsigned int m1, unsigned int n1,
150  Real alpha1, const Real* A)
151  {
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);
155 
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]);
157  }
158 };
159 
170 template <typename Real> void calc_f(Real* F, int n, Real t)
171 {
172  static ExpensiveMath expmath;
173  int i, m, k;
174  int m2;
175  Real t2;
176  Real num;
177  Real sum;
178  Real term1, term2;
179  static Real K = 1.0/M_2_SQRTPI;
180  Real et;
181 
182  if (t>20.0){ /* For big t's do upward recursion */
183  t2 = 2*t;
184  et = exp(-t);
185  t = sqrt(t);
186  F[0] = K*erf(t)/t;
187  for(m=0; m<=n-1; m++){
188  F[m+1] = ((2*m + 1)*F[m] - et)/(t2);
189  }
190  }
191  else { /* For smaller t's compute F with highest n using
192  asymptotic series (see I. Shavitt in
193  Methods in Computational Physics, ed. B. Alder eta l,
194  vol 2, 1963, page 8) */
195  et = exp(-t);
196  t2 = 2*t;
197  m2 = 2*n;
198  num = expmath.df[m2];
199  i=0;
200  sum = 1.0/(m2+1);
201  do{
202  i++;
203  num = num*t2;
204  term1 = num/expmath.df[m2+2*i+2];
205  sum += term1;
206  } while (term1 > EPS && i < expmath.MAXFAC);
207  F[n] = sum*et;
208  for(m=n-1;m>=0;m--){ /* And then do downward recursion */
209  F[m] = (t2*F[m+1] + et)/(2*m+1);
210  }
211  }
212 }
213 
214 namespace {
215  template <typename Int>
216  int parity(Int a) {
217  if (a%2 == 1) return -1;
218  return 1;
219  }
220 }
221 
222 namespace libint2 {
223  int min(int t1, unsigned int t2) { return t1 > (int)t2 ? (int)t2 : t1; }
224 };
225 
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,
244  int norm_flag)
245 {
246  static ExpensiveMath expmath;
247 
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]);
259 
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]);
272 
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;
277 
278  const LIBINT2_REF_REALTYPE four = 4;
279 
280  int u1,u2,v1,v2,w1,w2,tx,ty,tz,txmax,tymax,tzmax;
281  int i,j,k;
282  int lp,lq,mp,mq,np,nq;
283  int zeta;
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;
292 
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));
298 
299  if (norm_flag > 0) {
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);
304  }
305 
306  const int ltot = l1+l2+l3+l4+m1+m2+m3+m4+n1+n2+n3+n4;
307  F = new LIBINT2_REF_REALTYPE[ltot+1];
308  //calc_f<LIBINT2_REF_REALTYPE>(F,ltot,PQ2*gammapq);
310 
311 #define DEBUG 0
312 #if DEBUG
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;
317 #endif
318 
319  flp = new LIBINT2_REF_REALTYPE[l1+l2+1];
320  for(k=0;k<=l1+l2;k++) {
321  flp[k] = 0.0;
322  for(i=0;i<=libint2::min(k,l1);i++) {
323  j = k-i;
324  if (j > l2) continue;
325  tmp = expmath.bc[l1][i]*expmath.bc[l2][j];
326  if (l1 - i > 0)
327  tmp *= pow(PAx, l1 - i);
328  if (l2 - j > 0)
329  tmp *= pow(PBx, l2 - j);
330  flp[k] += tmp;
331  }
332  }
333  fmp = new LIBINT2_REF_REALTYPE[m1+m2+1];
334  for(k=0;k<=m1+m2;k++) {
335  fmp[k] = 0.0;
336  for(i=0;i<=libint2::min(k,m1);i++) {
337  j = k-i;
338  if (j > m2) continue;
339  tmp = expmath.bc[m1][i]*expmath.bc[m2][j];
340  if (m1 - i > 0)
341  tmp *= pow(PAy, m1 - i);
342  if (m2 - j > 0)
343  tmp *= pow(PBy, m2 - j);
344  fmp[k] += tmp;
345  }
346  }
347  fnp = new LIBINT2_REF_REALTYPE[n1+n2+1];
348  for(k=0;k<=n1+n2;k++) {
349  fnp[k] = 0.0;
350  for(i=0;i<=libint2::min(k,n1);i++) {
351  j = k-i;
352  if (j > n2) continue;
353  tmp = expmath.bc[n1][i]*expmath.bc[n2][j];
354  if (n1 - i > 0)
355  tmp *= pow(PAz, n1 - i);
356  if (n2 - j > 0)
357  tmp *= pow(PBz, n2 - j);
358  fnp[k] += tmp;
359  }
360  }
361  flq = new LIBINT2_REF_REALTYPE[l3+l4+1];
362  for(k=0;k<=l3+l4;k++) {
363  flq[k] = 0.0;
364  for(i=0;i<=libint2::min(k,l3);i++) {
365  j = k-i;
366  if (j > l4) continue;
367  tmp = expmath.bc[l3][i]*expmath.bc[l4][j];
368  if (l3 - i > 0)
369  tmp *= pow(QCx, l3 - i);
370  if (l4 - j > 0)
371  tmp *= pow(QDx, l4 - j);
372  flq[k] += tmp;
373  }
374  }
375  fmq = new LIBINT2_REF_REALTYPE[m3+m4+1];
376  for(k=0;k<=m3+m4;k++) {
377  fmq[k] = 0.0;
378  for(i=0;i<=libint2::min(k,m3);i++) {
379  j = k-i;
380  if (j > m4) continue;
381  tmp = expmath.bc[m3][i]*expmath.bc[m4][j];
382  if (m3 - i > 0)
383  tmp *= pow(QCy, m3 - i);
384  if (m4 - j > 0)
385  tmp *= pow(QDy, m4 - j);
386  fmq[k] += tmp;
387  }
388  }
389  fnq = new LIBINT2_REF_REALTYPE[n3+n4+1];
390  for(k=0;k<=n3+n4;k++) {
391  fnq[k] = 0.0;
392  for(i=0;i<=libint2:: min(k,n3);i++) {
393  j = k-i;
394  if (j > n4) continue;
395  tmp = expmath.bc[n3][i]*expmath.bc[n4][j];
396  if (n3 - i > 0)
397  tmp *= pow(QCz, n3 - i);
398  if (n4 - j > 0)
399  tmp *= pow(QDz, n4 - j);
400  fnq[k] += tmp;
401  }
402  }
403 
404  for (lp = 0; lp <= l1 + l2; lp++) {
405  for (lq = 0; lq <= l3 + l4; lq++) {
406  u1max = lp / 2;
407  u2max = lq / 2;
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++) {
418  v1max = mp / 2;
419  v2max = mq / 2;
420  for (v1 = 0; v1 <= v1max; v1++) {
421  for (v2 = 0; v2 <= v2max; v2++) {
422  Gy =
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++) {
433  w1max = np / 2;
434  w2max = nq / 2;
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
453  - tx - ty - tz;
454  result += Gx * Gy * Gz * F[zeta]
455  * parity(tx + ty + tz)
456  * pow(PQx,
457  lp + lq - 2 * u1 - 2 * u2 - 2 * tx)
458  * pow(PQy,
459  mp + mq - 2 * v1 - 2 * v2 - 2 * ty)
460  * pow(PQz,
461  np + nq - 2 * w1 - 2 * w2 - 2 * tz)
462  / (pow(
463  four,
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]);
473 #if DEBUG
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;
478 #endif
479  }
480  }
481  }
482  }
483  }
484  }
485  }
486  }
487  }
488  }
489  }
490  }
491  }
492  }
493  }
494 
495  delete[] F;
496  delete[] flp;
497  delete[] fmp;
498  delete[] fnp;
499  delete[] flq;
500  delete[] fmq;
501  delete[] fnq;
502 
503  return result*pfac;
504 }
505 
506 
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,
519  int norm_flag) {
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,
526  norm_flag);
527  else {
528 
529  struct nonzero_t {
530  bool operator()(unsigned int val) {
531  return val > 0;
532  }
533  };
534 
535  nonzero_t nonzero;
536  // handle one derivative at a time
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;
541 
542  unsigned int qn[4][3];
543  qn[0][0] = l1;
544  qn[0][1] = m1;
545  qn[0][2] = n1;
546  qn[1][0] = l2;
547  qn[1][1] = m2;
548  qn[1][2] = n2;
549  qn[2][0] = l3;
550  qn[2][1] = m3;
551  qn[2][2] = n3;
552  qn[3][0] = l4;
553  qn[3][1] = m4;
554  qn[3][2] = n4;
555 
556  LIBINT2_REF_REALTYPE alpha[4];
557  alpha[0] = alpha1;
558  alpha[1] = alpha2;
559  alpha[2] = alpha3;
560  alpha[3] = alpha4;
561 
562  ++qn[dc][dxyz];
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,
570  norm_flag) *
571  2.0 * alpha[dc];
572  --qn[dc][dxyz];
573 
574  if (qn[dc][dxyz] > 0) {
575  --qn[dc][dxyz];
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,
581  norm_flag) *
582  (qn[dc][dxyz] + 1);
583  ++qn[dc][dxyz];
584  }
585  return result;
586  }
587  assert(false); // unreachable;
588 }
589 
590 
591 
592 
593 #endif // header guard
Definition: eri.h:104
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