diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index 857ebc4..99b84f3 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -217,84 +217,27 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center)
   
 end function
  
-!double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
-! implicit none
-! include 'constants.F'
-!
-! real, intent(in)    :: a,b,gmma ! Exponents
-! real, intent(in)    :: xa,xb,xr ! Centers
-! integer, intent(in) :: i,j   ! Powers of xa and xb
-! integer :: ii, jj, kk, ll
-! real :: xp1,xp
-! real :: p1,p
-! double precision :: S(0:i+1,0:j+1)
-! double precision :: inv_p, di(max(i,j)), dj(j), c, c1
-!
-! ASSERT (a>0.)
-! ASSERT (b>0.)
-! ASSERT (i>=0)
-! ASSERT (j>=0)
-!
-! ! Gaussian product
-! call gaussian_product(a,xa,b,xb,c1,p1,xp1)
-! call gaussian_product(p1,xp1,gmma,xr,c,p,xp)
-! inv_p = 1.d0/p
-! S(0,0) = dsqrt(pi*inv_p)*c*c1
-! 
-! ! Obara-Saika recursion
-!
-! do ii=1,max(i,j)
-!  di(ii) = 0.5d0*inv_p*dble(ii)
-! enddo
-!
-! S(1,0) = (xp-xa) * S(0,0)
-! if (i>1) then
-!  do ii=1,i-1
-!   S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0)
-!  enddo
-! endif
-!
-! S(0,1) = (xp-xb) * S(0,0)
-! if (j>1) then
-!  do jj=1,j-1
-!   S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1)
-!  enddo
-! endif
-!
-! do jj=1,j
-!  S(1,jj) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1)
-!  do ii=2,i
-!   S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
-!  enddo
-! enddo
-!
-! ao_eplf_integral_primitive_oneD = S(i,j)
-!
-!end function
-
-double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
- implicit none
- include 'constants.F'
-
- real, intent(in)    :: a,b,gmma ! Exponents
- real, intent(in)    :: xa,xb,xr ! Centers
- integer, intent(in) :: i,j   ! Powers of xa and xb
- integer :: ii, jj, kk, ll
- real :: xp1,xp
- real :: p1,p
- double precision :: xpa, xpb
- double precision :: inv_p(2),S00, c
- double precision :: ObaraS
-
- ASSERT (a>0.)
- ASSERT (b>0.)
- ASSERT (i>=0)
- ASSERT (j>=0)
-
- ! Gaussian products
-!double precision :: t(2), xab(2), ab(2)
+ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
+  implicit none
+  include 'constants.F'
+ 
+  real, intent(in)    :: a,b,gmma ! Exponents
+  real, intent(in)    :: xa,xb,xr ! Centers
+  integer, intent(in) :: i,j   ! Powers of xa and xb
+  integer :: ii, jj, kk, ll
+  real :: xp1,xp
+  real :: p1,p
+  double precision :: S(0:i+1,0:j+1)
+  double precision :: inv_p(2), di(max(i,j)), dj(j), c
+ 
+  ASSERT (a>0.)
+  ASSERT (b>0.)
+  ASSERT (i>=0)
+  ASSERT (j>=0)
+ 
+  ! Gaussian product
+ ! Inlined Gaussian products (same as call gaussian_product)
  real :: t(2), xab(2), ab(2)
-!call gaussian_product(a,xa,b,xb,c1,p1,xp1)
  inv_p(1) = 1.d0/(a+b)
  p1 = a+b
  ab(1) = a*b
@@ -310,60 +253,136 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
  c = exp(- real(ab(1)*inv_p(1)*xab(1)**2 + &
              ab(2)*inv_p(2)*xab(2)**2) )
 
+! inv_p = 1.d0/p
+  S(0,0) = dsqrt(pi*inv_p(2))*c
+  
+  ! Obara-Saika recursion
+ 
+  do ii=1,max(i,j)
+   di(ii) = 0.5d0*inv_p(2)*dble(ii)
+  enddo
+ 
+  S(1,0) = (xp-xa) * S(0,0)
+  if (i>1) then
+   do ii=1,i-1
+    S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0)
+   enddo
+  endif
+ 
+  S(0,1) = (xp-xb) * S(0,0)
+  if (j>1) then
+   do jj=1,j-1
+    S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1)
+   enddo
+  endif
+ 
+  do jj=1,j
+   S(1,jj) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1)
+   do ii=2,i
+    S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
+   enddo
+  enddo
+ 
+  ao_eplf_integral_primitive_oneD = S(i,j)
+ 
+ end function
 
- xpa = xp-xa
- xpb = xp-xb
- S00 = dsqrt(pi*inv_p(2))*c
- ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p(2),S00)
-
-end function
-
-recursive double precision function ObaraS(i,j,xpa,xpb,inv_p,S00) result(res)
- implicit none
- integer, intent(in) :: i, j
- double precision, intent(in) :: xpa, xpb, inv_p
- double precision,intent(in) :: S00
-
- if (i == 0) then
-  if (j == 0) then
-   res = S00
-  else ! (j>0)
-   res = xpb*ObaraS(0,j-1,xpa,xpb,inv_p,S00)
-   if (j>1) then
-    res += 0.5d0*dble(j-1)*inv_p*ObaraS(0,j-2,xpa,xpb,inv_p,S00)
-   endif
-  endif ! (i==0).and.(j>0)
- else   ! (i>0)
-  if (j==0) then
-   res = xpa*ObaraS(i-1,0,xpa,xpb,inv_p,S00) 
-   if (i>1) then
-     res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,0,xpa,xpb,inv_p,S00) 
-   endif
-  else  ! (i>0).and.(j>0)
-   res = xpa * ObaraS(i-1,j,xpa,xpb,inv_p,S00) 
-   if (i>1) then
-    res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00) 
-   endif
-   res += 0.5d0*dble(j)*inv_p*ObaraS(i-1,j-1,xpa,xpb,inv_p,S00) 
-  endif  ! (i>0).and.(j>0)
- endif  ! (i>0)
-
-end function
+!double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
+! implicit none
+! include 'constants.F'
+!!
+! real, intent(in)    :: a,b,gmma ! Exponents
+! real, intent(in)    :: xa,xb,xr ! Centers
+! integer, intent(in) :: i,j   ! Powers of xa and xb
+! integer :: ii, jj, kk, ll
+! real :: xp1,xp
+! real :: p1,p
+! double precision :: xpa, xpb
+! double precision :: inv_p(2),S00, c
+! double precision :: ObaraS
+!!
+! ASSERT (a>0.)
+! ASSERT (b>0.)
+! ASSERT (i>=0)
+! ASSERT (j>=0)
+!!
+! ! Inlined Gaussian products (same as call gaussian_product)
+! real :: t(2), xab(2), ab(2)
+! inv_p(1) = 1.d0/(a+b)
+! p1 = a+b
+! ab(1) = a*b
+! inv_p(2) = 1.d0/(p1+gmma)
+! t(1) = (a*xa+b*xb)
+! xab(1) = xa-xb
+! xp1 = t(1)*inv_p(1)
+! p = p1+gmma
+! ab(2) = p1*gmma
+! t(2) = (p1*xp1+gmma*xr)
+! xab(2) = xp1-xr
+! xp = t(2)*inv_p(2)
+! c = exp(- real(ab(1)*inv_p(1)*xab(1)**2 + &
+!             ab(2)*inv_p(2)*xab(2)**2) )
+!!
+! xpa = xp-xa
+! xpb = xp-xb
+! S00 = sqrt(real(pi*inv_p(2)))*c
+! ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p(2),S00)
+!!
+!end function
+!!
+!recursive double precision function ObaraS(i,j,xpa,xpb,inv_p,S00) result(res)
+! implicit none
+! integer, intent(in) :: i, j
+! double precision, intent(in) :: xpa, xpb, inv_p
+! double precision,intent(in) :: S00
+!!
+! if (i == 0) then
+!  if (j == 0) then
+!   res = S00
+!  else ! (j>0)
+!   res = xpb*ObaraS(0,j-1,xpa,xpb,inv_p,S00)
+!   if (j>1) then
+!    res += 0.5d0*dble(j-1)*inv_p*ObaraS(0,j-2,xpa,xpb,inv_p,S00)
+!   endif
+!  endif ! (i==0).and.(j>0)
+! else   ! (i>0)
+!  if (j==0) then
+!   res = xpa*ObaraS(i-1,0,xpa,xpb,inv_p,S00) 
+!   if (i>1) then
+!     res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,0,xpa,xpb,inv_p,S00) 
+!   endif
+!  else  ! (i>0).and.(j>0)
+!   res = xpa * ObaraS(i-1,j,xpa,xpb,inv_p,S00) 
+!   if (i>1) then
+!    res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00) 
+!   endif
+!   res += 0.5d0*dble(j)*inv_p*ObaraS(i-1,j-1,xpa,xpb,inv_p,S00) 
+!  endif  ! (i>0).and.(j>0)
+! endif  ! (i>0)
+!!
+!end function
 
 double precision function ao_eplf_integral(i,j,gmma,center)
  implicit none
  integer, intent(in) :: i, j
  integer :: p,q,k
  double precision :: integral
+!DEC$ ATTRIBUTES FORCEINLINE
  double precision :: ao_eplf_integral_primitive_oneD
  real :: gmma, center(3)
+ double precision :: buffer(100)
+
 
  ASSERT(i>0)
  ASSERT(j>0)
  ASSERT(i<=ao_num)
  ASSERT(j<=ao_num)
 
- ao_eplf_integral = 0.
+ ao_eplf_integral = 0.d0
+ do p=1,ao_prim_num_max
+   buffer(p) = 0.d0
+ enddo
+
  do q=1,ao_prim_num(j)
   do p=1,ao_prim_num(i)
    integral =                              &
@@ -394,9 +413,13 @@ double precision function ao_eplf_integral(i,j,gmma,center)
     ao_power(j,3),                         &
     gmma,                                  &
     center(3))
-   ao_eplf_integral += integral*ao_coef(i,p)*ao_coef(j,q)
+!  ao_eplf_integral += integral*ao_coef(i,p)*ao_coef(j,q)
+   buffer(p) += integral*ao_coef(i,p)*ao_coef(j,q)
   enddo
  enddo
+ do p=1,ao_prim_num_max
+   ao_eplf_integral += buffer(p)
+ enddo
  
 end function