From 0e1c2f74a9e8eec073d01abb8f5f9affd06417c1 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@users.sourceforge.net>
Date: Tue, 5 Oct 2010 19:07:42 +0200
Subject: [PATCH] optimizations

---
 src/eplf_function.irp.f | 74 +++++++++++++++++++++++++----------------
 1 file changed, 46 insertions(+), 28 deletions(-)

diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index 6320a3f..67976d1 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -237,8 +237,8 @@ BEGIN_PROVIDER [ real, eplf_value_p ]
  if ( (aa > 0.d0).and.(ab > 0.d0) ) then
    aa = min(1.d0,aa)
    ab = min(1.d0,ab)
-   aa = -(dlog(aa)/eplf_gamma)
-   ab = -(dlog(ab)/eplf_gamma)
+   aa = -dlog(aa)/eplf_gamma
+   ab = -dlog(ab)/eplf_gamma
    aa = dsqrt(aa)
    ab = dsqrt(ab)
    eplf_value_p = (aa-ab)/(aa+ab+eps)
@@ -330,6 +330,7 @@ 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'
@@ -363,10 +364,15 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
  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) )
-
-  S(0,0) = dsqrt(pi*inv_p(2))*c
+ c = real(ab(1)*inv_p(1)*xab(1)**2 + &
+             ab(2)*inv_p(2)*xab(2)**2)
+ if ( c > 32.d0 ) then
+  ao_eplf_integral_primitive_oneD = 0.d0
+  return
+ endif
+ c = exp(-c) 
+ !S(0,0) = dsqrt(pi*inv_p(2))*c  
+  S(0,0) = 1.d0   ! Factor is applied at the end
   
   ! Obara-Saika recursion
  
@@ -397,7 +403,7 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
    enddo
   enddo
  
-  ao_eplf_integral_primitive_oneD = S(i,j)
+  ao_eplf_integral_primitive_oneD = dsqrt(pi*inv_p(2))*S(i,j)*c ! Application of the factor of S(0,0)
  
  end function
 
@@ -508,27 +514,32 @@ double precision function ao_eplf_integral(i,j,gmma,center)
     nucl_coord(ao_nucl(j),1),              &
     ao_power(j,1),                         &
     gmma,                                  &
-    center(1))   *                         &
-   ao_eplf_integral_primitive_oneD(        &
-    ao_expo(i,p),                          &
-    nucl_coord(ao_nucl(i),2),              &
-    ao_power(i,2),                         &
-    ao_expo(j,q),                          &
-    nucl_coord(ao_nucl(j),2),              &
-    ao_power(j,2),                         &
-    gmma,                                  &
-    center(2))   *                         &
-   ao_eplf_integral_primitive_oneD(        &
-    ao_expo(i,p),                          &
-    nucl_coord(ao_nucl(i),3),              &
-    ao_power(i,3),                         &
-    ao_expo(j,q),                          &
-    nucl_coord(ao_nucl(j),3),              &
-    ao_power(j,3),                         &
-    gmma,                                  &
-    center(3))
-!  ao_eplf_integral += integral*ao_coef(i,p)*ao_coef(j,q)
-   buffer(p) += integral*ao_coef(i,p)*ao_coef(j,q)
+    center(1))
+   if (integral /= 0.d0) then
+     integral *=                             &
+      ao_eplf_integral_primitive_oneD(       &
+      ao_expo(i,p),                          &
+      nucl_coord(ao_nucl(i),2),              &
+      ao_power(i,2),                         &
+      ao_expo(j,q),                          &
+      nucl_coord(ao_nucl(j),2),              &
+      ao_power(j,2),                         &
+      gmma,                                  &
+      center(2))
+     if (integral /= 0.d0) then
+      integral *=                             &
+       ao_eplf_integral_primitive_oneD(       &
+       ao_expo(i,p),                          &
+       nucl_coord(ao_nucl(i),3),              &
+       ao_power(i,3),                         &
+       ao_expo(j,q),                          &
+       nucl_coord(ao_nucl(j),3),              &
+       ao_power(j,3),                         &
+       gmma,                                  &
+       center(3))
+      buffer(p) += integral*ao_coef(i,p)*ao_coef(j,q)
+     endif
+   endif
   enddo
  enddo
  do p=1,ao_prim_num_max
@@ -555,3 +566,10 @@ double precision function mo_eplf_integral(i,j)
  
 end function
 
+
+
+
+
+
+
+