From 5de7b720283db5836f3d6023141669bd250a5c6e Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 20 May 2009 17:57:35 +0200
Subject: [PATCH] Added recursive fortran func

---
 Makefile   |   5 +--
 eplf.irp.f | 101 ++++++++++++++++++++++++++++++++++++++++++++++-------
 2 files changed, 92 insertions(+), 14 deletions(-)

diff --git a/Makefile b/Makefile
index bf3167b..87e8ff6 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,7 @@
 IRPF90 = irpf90 -DMPI#-a -d 
-FC     = mpif90 
-FCFLAGS= -O3 -xT
+IRPF90 = irpf90 -DMPI#-a -d 
+FC     = mpif90 -xT
+FCFLAGS= -O3 -g
 
 #FC     = gfortran -g -ffree-line-length-none
 #FCFLAGS=
diff --git a/eplf.irp.f b/eplf.irp.f
index 4e98347..01981bf 100644
--- a/eplf.irp.f
+++ b/eplf.irp.f
@@ -33,13 +33,15 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
   do j=i,mo_num
     mo_eplf_integral_matrix(j,i) = 0.
   enddo
+ enddo
 
-  do k=1,ao_num
+  do i=1,mo_num
+   do k=1,ao_num
     if (mo_coef(k,i) /= 0.) then
      do j=i,mo_num
        do l=1,ao_num
         mo_eplf_integral_matrix(j,i) = mo_eplf_integral_matrix(j,i) + &
-          mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l)
+          mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(l,k)
        enddo
      enddo
     endif
@@ -203,7 +205,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)
+double precision function ao_eplf_integral_primitive_oneD2(a,xa,i,b,xb,j,gmma,xr)
  implicit none
  include 'constants.F'
 
@@ -213,7 +215,7 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
  integer :: ii, jj, kk, ll
  real :: xp1,xp
  real :: p1,p
- double precision :: S(0:i+1,0:j)
+ double precision :: S(0:i+1,0:j+1)
  double precision :: inv_p, di(max(i,j)), dj(j), c, c1
 
  ASSERT (a>0.)
@@ -229,23 +231,18 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
  
  ! Obara-Saika recursion
 
- if (i>0) then
-   S(1,0) = (xp-xa) * S(0,0)
- endif
- if (j>0) then
-   S(0,1) = (xp-xb) * S(0,0)
- endif
-
  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)
@@ -259,7 +256,69 @@ 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_oneD2 = 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 :: S00, xpa, xpb
+ double precision :: inv_p,c,c1
+ double precision :: ObaraS
+
+ 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
+ S00 = dsqrt(pi*inv_p)*c*c1
+ xpa = xp-xa
+ xpb = xp-xb
+ ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p,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 = 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 = 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 = res + 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00) 
+   endif
+   res = 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
 
@@ -313,3 +372,21 @@ double precision function ao_eplf_integral(i,j,gmma,center)
  
 end function
 
+ double precision function mo_eplf_integral(i,j)
+  implicit none
+  integer :: i, j, k, l
+  PROVIDE ao_eplf_integral_matrix
+  PROVIDE mo_coef
+  mo_eplf_integral = 0.d0
+ 
+   do k=1,ao_num
+     if (mo_coef(k,i) /= 0.) then
+       do l=1,ao_num
+        mo_eplf_integral = mo_eplf_integral + &
+          mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l)
+       enddo
+     endif
+   enddo
+ 
+ end function
+