BEGIN_PROVIDER  [ real, eplf_gamma ]
 implicit none
 BEGIN_DOC  
! Value of the gaussian for the EPLF
 END_DOC
 include 'constants.F'
 real            :: eps, N
! N = 0.1
! eps = -real(dlog(tiny(1.d0)))
! eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps
 N=0.001
 eps = 50.
 eplf_gamma = eps**2 *0.5d0*(4./9.*pi*density_value_p * N**(-1./3.))**(2./3.)
!eplf_gamma = 1.e10
!eplf_gamma = 1.e5
END_PROVIDER

BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
 implicit none
 BEGIN_DOC  
! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
 END_DOC
 integer :: i, j
 double precision :: ao_eplf_integral
 do i=1,ao_num
  do j=i,ao_num
    ao_eplf_integral_matrix(j,i) = ao_eplf_integral(j,i,dble(eplf_gamma),point)
    ao_eplf_integral_matrix(i,j) = ao_eplf_integral_matrix(j,i)
  enddo
 enddo
END_PROVIDER

BEGIN_PROVIDER [ integer, ao_eplf_integral_matrix_non_zero_idx, (0:ao_num,ao_num) ]
 implicit none
 BEGIN_DOC  
! Non zero indices of the ao_eplf_integral_matrix
 END_DOC
 integer :: i, j
 integer :: idx
 do j=1,ao_num
  idx = 0
  do i=1,ao_num
   if (abs(ao_eplf_integral_matrix(i,j)) > 1.e-16) then
    idx +=1
    ao_eplf_integral_matrix_non_zero_idx(idx,j) = i
   endif
  enddo
  ao_eplf_integral_matrix_non_zero_idx(0,j) = idx
 enddo
END_PROVIDER

BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
 implicit none
 BEGIN_DOC  
! Array of all the <phi_i phi_j | exp(-gamma r^2)> for EPLF
 END_DOC
 integer :: i, j, k, l, kmo, kao
 double precision :: t, moki
 !DEC$ VECTOR ALIGNED
 mo_eplf_integral_matrix = 0.d0

 do k=1,ao_num
  do kmo=1,mo_coef_transp_non_zero_idx(0,k)
   i = mo_coef_transp_non_zero_idx(kmo,k)
   moki = mo_coef_transp(i,k)
   do kao = 1,ao_eplf_integral_matrix_non_zero_idx(0,k)
     l = ao_eplf_integral_matrix_non_zero_idx(kao,k)
     t = moki*ao_eplf_integral_matrix(l,k)
     do j=i,mo_num
       mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l)
     enddo
   enddo
  enddo

 enddo

!do j=1,mo_num
! do i=1,j-1
!   mo_eplf_integral_matrix(i,j) = mo_eplf_integral_matrix(j,i)
! enddo
!enddo

 integer :: momax
 momax = mo_num - mod(mo_num,4)
 do j=1,momax,4
  do i=1,j-1
    mo_eplf_integral_matrix(i,j  ) = mo_eplf_integral_matrix(j  ,i)
    mo_eplf_integral_matrix(i,j+1) = mo_eplf_integral_matrix(j+1,i)
    mo_eplf_integral_matrix(i,j+2) = mo_eplf_integral_matrix(j+2,i)
    mo_eplf_integral_matrix(i,j+3) = mo_eplf_integral_matrix(j+3,i)
  enddo
  mo_eplf_integral_matrix(j  ,j+1) = mo_eplf_integral_matrix(j+1,j)
  mo_eplf_integral_matrix(j  ,j+2) = mo_eplf_integral_matrix(j+2,j)
  mo_eplf_integral_matrix(j  ,j+3) = mo_eplf_integral_matrix(j+3,j)
  mo_eplf_integral_matrix(j+1,j+2) = mo_eplf_integral_matrix(j+2,j+1)
  mo_eplf_integral_matrix(j+1,j+3) = mo_eplf_integral_matrix(j+3,j+1)
  mo_eplf_integral_matrix(j+2,j+3) = mo_eplf_integral_matrix(j+3,j+2)
 enddo

 do j=momax+1,mo_num
  do i=1,j-1
    mo_eplf_integral_matrix(i,j) = mo_eplf_integral_matrix(j,i)
  enddo
 enddo

 
END_PROVIDER


 BEGIN_PROVIDER [ double precision, eplf_up_up ]
&BEGIN_PROVIDER [ double precision, eplf_up_dn ]
 implicit none
 BEGIN_DOC  
! Value of the d_upup and d_updn quantities needed for EPLF
 END_DOC

 integer :: i, j

 eplf_up_up = 0.d0
 eplf_up_dn = 0.d0

 PROVIDE mo_coef_transp

 do i=1,mo_closed_num
  do j=1,mo_closed_num
    double precision :: temp
    temp = mo_value_prod_p(i,i)*mo_eplf_integral_matrix(j,j)
    eplf_up_up += temp - mo_value_prod_p(j,i)*mo_eplf_integral_matrix(j,i) 
    eplf_up_dn += temp
  enddo
 enddo
 eplf_up_up *= 2.d0
 eplf_up_dn *= 2.d0 

 integer :: k,l,m

 do m=1,two_e_density_num
  i=two_e_density_indice(1,m)
  j=two_e_density_indice(2,m)
  k=two_e_density_indice(3,m)
  l=two_e_density_indice(4,m)
  temp = mo_value_prod_p(i,j)*mo_eplf_integral_matrix(k,l)
  eplf_up_up += two_e_density_value(1,m)*temp
  eplf_up_dn += two_e_density_value(2,m)*temp
 enddo

END_PROVIDER


BEGIN_PROVIDER [ real, eplf_value_p ]
 implicit none
 BEGIN_DOC  
! Value of the EPLF at the current point.
 END_DOC
 double precision :: aa, ab
 double precision, parameter :: eps = tiny(1.d0)

 aa = eplf_up_up
 ab = eplf_up_dn
 if ( (aa > 0.d0).and.(ab > 0.d0) ) then
   aa = min(1.d0,aa)
   ab = min(1.d0,ab)
   aa = -dlog(aa)
   aa = dsqrt(aa)
   ab = -dlog(ab)
   ab = dsqrt(ab)
   eplf_value_p = (aa-ab)/(aa+ab+eps)
 else
   eplf_value_p = 0.d0
 endif

END_PROVIDER


!double precision function ao_eplf_integral_primitive_oneD_numeric(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,parameter :: Npoints=10000
! real :: x, xmin, xmax, dx
!
! ASSERT (a>0.)
! ASSERT (b>0.)
! ASSERT (i>=0)
! ASSERT (j>=0)
!
! xmin = min(xa,xb)
! xmax = max(xa,xb)
! xmin = min(xmin,xr) - 10.
! xmax = max(xmax,xr) + 10.
! dx = (xmax-xmin)/real(Npoints)
!
! real :: dtemp
! dtemp = 0.
! x = xmin
! integer :: k
! do k=1,Npoints
!  dtemp += &
!   (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2))
!  x = x+dx
! enddo
! ao_eplf_integral_primitive_oneD_numeric = dtemp*dx
!
!end function
!
!double precision function ao_eplf_integral_numeric(i,j,gmma,center)
! implicit none
! integer, intent(in) :: i, j
! integer :: p,q,k
! double precision :: integral
! double precision :: ao_eplf_integral_primitive_oneD_numeric
! real :: gmma, center(3), c
!
! 
!  ao_eplf_integral_numeric = 0.d0
!  do q=1,ao_prim_num(j)
!   do p=1,ao_prim_num(i)
!    c = ao_coef(p,i)*ao_coef(q,j)
!    integral =                              &
!    ao_eplf_integral_primitive_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),                         &
!     gmma,                                  &
!     center(1))   *                         &
!    ao_eplf_integral_primitive_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),                         &
!     gmma,                                  &
!     center(2))   *                         &
!    ao_eplf_integral_primitive_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),                         &
!     gmma,                                  &
!     center(3))
!    ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral
!   enddo
!  enddo
!  
!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 ! Exponents
  double precision , intent(in) :: gmma  ! eplf_gamma 
  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.)
  ASSERT (i>=0)
  ASSERT (j>=0)
 
  double precision :: alpha, beta, delta, c, beta_inv
  double precision :: db_inv

  beta     = a       + b       + gmma
  beta_inv = 1.d0/beta
  delta    = a*xa    + b*xb    + gmma*xr
  db_inv = delta*beta_inv
  alpha    = a*xa*xa + b*xb*xb + gmma*xr*xr
  c = alpha-delta*db_inv
  
  if ( c > 20.d0 ) then   ! Cut-off on exp(-32)
   ao_eplf_integral_primitive_oneD = 0.d0
   return
  endif

  double precision :: u,v,w
  c = exp(-c) 
  w = sqrt(beta_inv)
  u=db_inv-xa
  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
   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
 
end function


double precision function ao_eplf_integral(i,j,gmma,center)
 implicit none
 integer, intent(in) :: i, j
 real, intent(in)    :: center(3)
 double precision, intent(in) :: gmma
 integer :: p,q,k
 double precision :: integral
 double precision :: ao_eplf_integral_primitive_oneD


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

 ao_eplf_integral = 0.d0

 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),                          &
        nucl_coord(ao_nucl(i),1),              &
        ao_power(i,1),                         &
        ao_expo(q,j),                          &
        nucl_coord(ao_nucl(j),1),              &
        ao_power(j,1),                         &
        gmma,                                  &
        center(1))
       if (dabs(integral) > 1.d-15)  then
         integral *=                             &
          ao_eplf_integral_primitive_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),                         &
          gmma,                                  &
          center(2))
         if (dabs(integral) > 1.d-15)  then
          integral *=                             &
           ao_eplf_integral_primitive_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),                         &
           gmma,                                  &
           center(3))
          ao_eplf_integral += integral*ao_coef(p,i)*ao_coef(q,j)
         endif
       endif
    enddo
 enddo
 
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_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(l,k)
       enddo
     endif
   enddo
 
end function