From 40ad59efa8283949c2b767432d171f140ff87b39 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@users.sourceforge.net>
Date: Thu, 31 May 2012 23:26:30 +0200
Subject: [PATCH] Optimized eplf integrals

---
 src/Util.irp.f          | 24 ++++++++++++----
 src/eplf_function.irp.f | 62 +++++++++++++++++++++++++++++++----------
 2 files changed, 66 insertions(+), 20 deletions(-)

diff --git a/src/Util.irp.f b/src/Util.irp.f
index 25076c5..2b5c884 100644
--- a/src/Util.irp.f
+++ b/src/Util.irp.f
@@ -23,7 +23,18 @@ double precision function binom(k,n)
  !DEC$ ATTRIBUTES ALIGN : 32 :: memo
  integer, save :: ifirst
  double precision :: fact
- if (ifirst == 0) then
+ if (ifirst /= 0) then
+
+  if (n<20) then
+    binom = memo(k,n)
+    return
+  else
+    binom=fact(n)/(fact(k)*fact(n-k))
+    return
+  endif
+
+ else
+
   ifirst = 1
   integer :: i,j
   do j=0,20
@@ -31,12 +42,13 @@ double precision function binom(k,n)
     memo(i,j) = fact(j)/(fact(i)*fact(j-i))
    enddo
   enddo
- endif
 
- if (n<20) then
-   binom = memo(k,n)
- else
-   binom=fact(n)/(fact(k)*fact(n-k))
+  if (n<20) then
+    binom = memo(k,n)
+    return
+  else
+    binom=fact(n)/(fact(k)*fact(n-k))
+  endif
  endif
 end
 
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index 9590890..c8125ef 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -243,6 +243,16 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
   real, intent(in)    :: xa,xb,xr ! Centers
   integer, intent(in) :: i,j   ! Powers of xa and xb
   integer :: ii, jj, kk
+
+  integer, save :: ifirst=0
+  double precision, save :: f2(32)
+
+  if (ifirst == 0) then
+    ifirst = 1
+    do ii=1,32
+      f2(ii) = 2.d0**(-ii/2) * fact2(ii-1)
+    enddo
+  endif
  
   ASSERT (a>0.)
   ASSERT (b>0.)
@@ -271,16 +281,43 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
   v=db_inv-xb
 
   double precision :: fact2, binom, f, ac
+  double precision :: ui(-1:i), vj(-1:j), wij(-1:i+j+3)
+  !DEC$ ATTRIBUTES ALIGN : 32 :: ui, vj, wij
+
+  ui(i) = 1.d0
+  ui(i-1) = u
+  do ii=i-2,0,-1
+   ui(ii) = ui(ii+1) * u
+  enddo
+
+  vj(j) = 1.d0
+  vj(j-1) = v
+  do jj=j-2,0,-1
+   vj(jj) = vj(jj+1) * v
+  enddo
+
+  wij(0) = 1.d0
+  wij(1) = w 
+  wij(2) = w*w
+  do ii=3,i+j
+   wij(ii) = wij(ii-1) * w 
+  enddo
+  
+  do ii=2,i+j
+   wij(ii) *= f2(ii) 
+  enddo
+  
   ac = 0.d0
   integer :: istart(0:1)
   istart(0) = 0
   istart(1) = 1
   do ii=0,i
-   do jj=istart(mod(ii,2)),j,2
-     kk= -(ii+jj)/2
-     f=binom(ii,i)*binom(jj,j)*fact2(ii+jj-1)
-     ac += u**(i-ii) * v**(j-jj) * w**(ii+jj) * f * 2.d0**kk
-   enddo
+   if (istart(mod(ii,2))<=j) then
+     f=ui(ii)*binom(ii,i)
+     do jj=istart(mod(ii,2)),j,2
+       ac += vj(jj) * wij(ii+jj) * f *binom(jj,j)
+     enddo
+   endif
   enddo
   ao_eplf_integral_primitive_oneD = dsqrt_pi*w*c*ac
  
@@ -305,12 +342,12 @@ double precision function ao_eplf_integral(i,j,gmma,center)
  ASSERT(j<=ao_num)
 
  ao_eplf_integral = 0.d0
- do p=1,ao_prim_num(i)
-   buffer(p) = 0.d0
- enddo
 
  do q=1,ao_prim_num(j)
     do p=1,ao_prim_num(i)
+       if (abs(ao_coef(q,j)*ao_coef(p,i)) < 1.e-9) then
+         cycle
+       endif
        integral =                              &
        ao_eplf_integral_primitive_oneD(        &
         ao_expo(p,i),                          &
@@ -321,7 +358,7 @@ double precision function ao_eplf_integral(i,j,gmma,center)
         ao_power(j,1),                         &
         gmma,                                  &
         center(1))
-       if (integral /= 0.d0) then
+       if (dabs(integral) > 1.d-15)  then
          integral *=                             &
           ao_eplf_integral_primitive_oneD(       &
           ao_expo(p,i),                          &
@@ -332,7 +369,7 @@ double precision function ao_eplf_integral(i,j,gmma,center)
           ao_power(j,2),                         &
           gmma,                                  &
           center(2))
-         if (integral /= 0.d0) then
+         if (dabs(integral) > 1.d-15)  then
           integral *=                             &
            ao_eplf_integral_primitive_oneD(       &
            ao_expo(p,i),                          &
@@ -343,14 +380,11 @@ double precision function ao_eplf_integral(i,j,gmma,center)
            ao_power(j,3),                         &
            gmma,                                  &
            center(3))
-          buffer(p) += integral*ao_coef(p,i)*ao_coef(q,j)
+          ao_eplf_integral += integral*ao_coef(p,i)*ao_coef(q,j)
          endif
        endif
     enddo
  enddo
- do p=1,ao_prim_num(i)
-   ao_eplf_integral += buffer(p)
- enddo
  
 end function