BEGIN_PROVIDER [ double precision, ao_overlap_matrix, (ao_num,ao_num) ]
 implicit none
 BEGIN_DOC  
! Overlap matrix between the Atomic Orbitals
 END_DOC

 integer :: i, j
 double precision :: ao_overlap
 do j=1,ao_num
  do i=1,j
   ao_overlap_matrix(i,j) = ao_overlap(i,j)
   ao_overlap_matrix(j,i) = ao_overlap_matrix(i,j)
  enddo
 enddo
END_PROVIDER


double precision function primitive_overlap_oneD_numeric(a,xa,i,b,xb,j)
 implicit none
 include 'constants.F'

 real, intent(in)    :: a,b   ! Exponents
 real, intent(in)    :: xa,xb ! Centers
 integer, intent(in) :: i,j   ! Powers of xa and xb
 integer,parameter :: Npoints=10000
 real :: x, xmin, xmax, dx

 ASSERT (a>0.)
 ASSERT (b>0.)
 ASSERT (i>=0)
 ASSERT (j>=0)

 xmin = min(xa,xb) - 10.
 xmax = max(xa,xb) + 10.
 dx = (xmax-xmin)/real(Npoints)

 real :: dtemp
 dtemp = 0.
 x = xmin
 integer :: k
 do k=1,Npoints
  dtemp = dtemp + &
   (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2))
  x = x+dx
 enddo
 primitive_overlap_oneD_numeric = dtemp*dx

end function

double precision function ao_overlap_numeric(i,j)
 implicit none
 integer, intent(in) :: i, j
 integer :: p,q,k
 double precision :: integral(ao_prim_num_max,ao_prim_num_max)

 double precision :: primitive_overlap_oneD_numeric

 ASSERT(i>0)
 ASSERT(j>0)
 ASSERT(i<=ao_num)
 ASSERT(j<=ao_num)

 do q=1,ao_prim_num(j)
  do p=1,ao_prim_num(i)
   integral(p,q) =                         &
   primitive_overlap_oneD_numeric (        &
    ao_expo(p,i),                          &
    nucl_coord(ao_nucl(i),1),              &
    ao_power(i,1),                         &
    ao_expo(q,j),                          &
    nucl_coord(ao_nucl(j),1),              &
    ao_power(j,1) ) *                      &
   primitive_overlap_oneD_numeric (        &
    ao_expo(p,i),                          &
    nucl_coord(ao_nucl(i),2),              &
    ao_power(i,2),                         &
    ao_expo(q,j),                          &
    nucl_coord(ao_nucl(j),2),              &
    ao_power(j,2) ) *                      &
   primitive_overlap_oneD_numeric (        &
    ao_expo(p,i),                          &
    nucl_coord(ao_nucl(i),3),              &
    ao_power(i,3),                         &
    ao_expo(q,j),                          &
    nucl_coord(ao_nucl(j),3),              &
    ao_power(j,3) ) 
  enddo
 enddo

 do q=1,ao_prim_num(j)
  do p=1,ao_prim_num(i)
   integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j)
  enddo
 enddo

 ao_overlap_numeric = 0.
 do q=1,ao_prim_num(j)
  do p=1,ao_prim_num(i)
   ao_overlap_numeric = ao_overlap_numeric + integral(p,q)
  enddo
 enddo
 
end function

subroutine gaussian_product(a,xa,b,xb,k,p,xp)
 implicit none
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}

 real, intent(in) :: a,b   ! Exponents
 real, intent(in) :: xa,xb ! Centers
 real, intent(out) :: p    ! New exponent
 real, intent(out) :: xp   ! New center
 double precision, intent(out) :: k  ! Constant

 double precision:: p_inv

 ASSERT (a>0.)
 ASSERT (b>0.)

 double precision :: t, xab, ab

 p_inv = 1.d0/(a+b)
 p = a+b
 ab = a*b
 t = (a*xa+b*xb)
 xab = xa-xb
 xp = t*p_inv
 k = dexp(-ab*p_inv*xab**2)

end subroutine

double precision function primitive_overlap_oneD(a,xa,i,b,xb,j)
 implicit none
 include 'constants.F'

 real, intent(in)    :: a,b   ! Exponents
 real, intent(in)    :: xa,xb ! Centers
 integer, intent(in) :: i,j   ! Powers of xa and xb
 integer :: ii, jj, kk, ll
 real :: xp
 real :: p
 double precision :: S(0:i+1,0:j)
 double precision :: inv_p, di(max(i,j)), dj(j), c

 ASSERT (a>0.)
 ASSERT (b>0.)
 ASSERT (i>=0)
 ASSERT (j>=0)

 ! Gaussian product
 p = a+b
 xp = (a*xa+b*xb)
 inv_p = 1.d0/p
 xp = xp*inv_p

 c = a*b*inv_p*(xa-xb)**2
 if ( c > 32.d0 ) then   ! Cut-off on exp(-32)
  primitive_overlap_oneD = 0.d0
  return
 endif

 c =  dexp(-c)

 S(0,0) = 1.d0

 ! Obara-Saika recursion

 if (i>0) then
   S(1,0) = xp-xa
 endif
 if (j>0) then
   S(0,1) = xp-xb
 endif

 do ii=1,max(i,j)
  di(ii) = 0.5d0*inv_p*dble(ii)
 enddo

 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

 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

 primitive_overlap_oneD = S(i,j)*c* dsqrt(pi*inv_p)

end function

double precision function ao_overlap(i,j)
 implicit none
 integer, intent(in) :: i, j
 integer :: p,q,k
 double precision :: integral(ao_prim_num_max,ao_prim_num_max)

 double precision :: primitive_overlap_oneD

 ASSERT(i>0)
 ASSERT(j>0)
 ASSERT(i<=ao_num)
 ASSERT(j<=ao_num)

 do q=1,ao_prim_num(j)
  do p=1,ao_prim_num(i)
   integral(p,q) =                         &
   primitive_overlap_oneD (                &
    ao_expo(p,i),                          &
    nucl_coord(ao_nucl(i),1),              &
    ao_power(i,1),                         &
    ao_expo(q,j),                          &
    nucl_coord(ao_nucl(j),1),              &
    ao_power(j,1) ) *                      &
   primitive_overlap_oneD (                &
    ao_expo(p,i),                          &
    nucl_coord(ao_nucl(i),2),              &
    ao_power(i,2),                         &
    ao_expo(q,j),                          &
    nucl_coord(ao_nucl(j),2),              &
    ao_power(j,2) ) *                      &
   primitive_overlap_oneD (                &
    ao_expo(p,i),                          &
    nucl_coord(ao_nucl(i),3),              &
    ao_power(i,3),                         &
    ao_expo(q,j),                          &
    nucl_coord(ao_nucl(j),3),              &
    ao_power(j,3) ) 
  enddo
 enddo

 do q=1,ao_prim_num(j)
  do p=1,ao_prim_num(i)
   integral(p,q) *= ao_coef(p,i)*ao_coef(q,j)
  enddo
 enddo

 ao_overlap = 0.d0
 do q=1,ao_prim_num(j)
  do p=1,ao_prim_num(i)
   ao_overlap += integral(p,q)
  enddo
 enddo
 
end function