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
!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 [ 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
 double precision :: t
 PROVIDE ao_eplf_integral_matrix
 PROVIDE mo_coef
 do i=1,mo_num
  do j=i,mo_num
    mo_eplf_integral_matrix(j,i) = 0.d0
  enddo


  do k=1,ao_num
   if (mo_coef(k,i) /= 0.) then
    do l=1,ao_num
     if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then
      t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k)
      do j=i,mo_num
       mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l)
      enddo
     endif
    enddo
   endif

  enddo

 enddo

 do i=1,mo_num
  do j=i+1,mo_num
    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
    eplf_up_up += mo_value_p(i)* (  &
      mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
      mo_value_p(j)*mo_eplf_integral_matrix(j,i) )
    eplf_up_dn += mo_value_p(i)*mo_value_p(i)* &
        mo_eplf_integral_matrix(j,j)
  enddo
 enddo
 eplf_up_up *= 2.d0
 eplf_up_dn *= 2.d0 

 integer :: k,l,m,n,p,p2
 integer :: ik,il,jk,jl
 double precision :: phase,dtemp(2)
 integer :: exc(4)

 PROVIDE det
 PROVIDE elec_num_2

 do k=1,det_num
  do l=1,det_num

   exc(1) = abs(det_exc(k,l,1))
   exc(2) = abs(det_exc(k,l,2))
   exc(3) = exc(1)+exc(2)
   exc(4) = det_exc(k,l,1)*det_exc(k,l,2)
   if (exc(4) /= 0) then
    exc(4) = exc(4)/abs(exc(4))
   else
    exc(4) = 1
   endif

   dtemp(1) = 0.d0
   dtemp(2) = 0.d0
   do p=1,2
     p2 = 1+mod(p,2)
     if ( exc(3) == 0 ) then
      ! Closed-open shell interactions
      do j=1,mo_closed_num
       do n=1,elec_num_2(p)-mo_closed_num
        ik = det(n,k,p)
        il = det(n,l,p)
        dtemp(1) += mo_value_p(ik)* (  &
           mo_value_p(il)*mo_eplf_integral_matrix(j,j) - &
           mo_value_p(j)*mo_eplf_integral_matrix(j,il) )
        dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(j,j) 
       enddo
      enddo

      !- Open-closed shell interactions
      do m=1,elec_num_2(p)-mo_closed_num
       jk = det(m,k,p)
       jl = det(m,l,p)
       do i=1,mo_closed_num
         dtemp(1) += mo_value_p(i)* (  &
           mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) - &
           mo_value_p(jl)*mo_eplf_integral_matrix(jk,i) )
         dtemp(2) += mo_value_p(i)*mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) 
       enddo
      enddo

      !- Open-open shell interactions
      do m=1,elec_num_2(p)-mo_closed_num
       jk = det(m,k,p)
       jl = det(m,l,p)
       do n=1,elec_num_2(p)-mo_closed_num
         ik = det(n,k,p)
         il = det(n,l,p)
         dtemp(1) += mo_value_p(ik)* (  &
           mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
           mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
       enddo
       do n=1,elec_num_2(p2)-mo_closed_num
         ik = det(n,k,p2)
         il = det(n,l,p2)
         dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) 
       enddo
      enddo
 
     else if ( (exc(3) == 1).and.(exc(p) == 1) ) then

     ! Sum over only the sigma-sigma interactions involving the excitation
      call get_single_excitation(k,l,ik,il,p)

      !- Open-closed shell interactions
      do j=1,mo_closed_num
       dtemp(1) += mo_value_p(ik)* (  &
           mo_value_p(il)*mo_eplf_integral_matrix(j,j) - &
           mo_value_p(j)*mo_eplf_integral_matrix(j,il) )
       dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(j,j)
      enddo

      !- Closed-open shell interactions
      do i=1,mo_closed_num
       dtemp(1) += mo_value_p(i)* (  &
           mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) - &
           mo_value_p(jl)*mo_eplf_integral_matrix(jk,i) )
       dtemp(2) += mo_value_p(i)*mo_value_p(i)*mo_eplf_integral_matrix(jk,jl)
      enddo

      !- Open-open shell interactions
      do m=1,elec_num_2(p)-mo_closed_num
       jk = det(m,k,p)
       jl = det(m,l,p)
       dtemp(1) += mo_value_p(ik)* (  &
           mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
           mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
      enddo
      do m=1,elec_num_2(p2)-mo_closed_num
       jk = det(m,k,p2)
       jl = det(m,l,p2)
       dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) 
      enddo

     else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
    
     ! Consider only the double excitations of same-spin electrons
      call get_double_excitation(k,l,ik,il,jk,jl,p)
    
      dtemp(1) += mo_value_p(ik)* (  &
             mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
             mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
    
     else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
    
     ! Consider only the double excitations of opposite-spin electrons
     call get_single_excitation(k,l,ik,il,p)
     call get_single_excitation(k,l,jk,jl,p2)
    
      dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl)
    
     endif
   enddo

   phase = dble(exc(4))
   eplf_up_up += phase * det_coef(k)*det_coef(l) * dtemp(1) 
   eplf_up_dn += phase * det_coef(k)*det_coef(l) * dtemp(2)

  enddo
 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)
   ab = -dlog(ab)
   aa = dsqrt(aa)
   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(i,p)*ao_coef(j,q)
!    integral =                              &
!    ao_eplf_integral_primitive_oneD_numeric(   &
!     ao_expo(i,p),                          &
!     nucl_coord(ao_nucl(i),1),              &
!     ao_power(i,1),                         &
!     ao_expo(j,q),                          &
!     nucl_coord(ao_nucl(j),1),              &
!     ao_power(j,1),                         &
!     gmma,                                  &
!     center(1))   *                         &
!    ao_eplf_integral_primitive_oneD_numeric(   &
!     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_numeric(   &
!     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_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, 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
  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 = dble(ab(1)*inv_p(1)*xab(1)**2 + &
              ab(2)*inv_p(2)*xab(2)**2)
  
! double precision, save :: c_accu(2)
! c_accu(1) += c
! c_accu(2) += 1.d0
! print *, c_accu(1)/c_accu(2)

  if ( c > 32.d0 ) then   ! Cut-off on exp(-32)
   ao_eplf_integral_primitive_oneD = 0.d0
   return
  endif

  c = exp(-c) 
  
  ! Obara-Saika recursion
  S(0,0) = 1.d0
 
  do ii=1,max(i,j)
   di(ii) = 0.5d0*inv_p(2)*dble(ii)
  enddo
 
  xab(1) = xp-xa
  xab(2) = xp-xb
  S(1,0) = xab(1) ! * S(0,0)
  if (i>1) then
   do ii=1,i-1
    S(ii+1,0) = xab(1) * S(ii,0) + di(ii)*S(ii-1,0)
   enddo
  endif
 
  S(0,1) = xab(2) ! * S(0,0)
  if (j>1) then
   do jj=1,j-1
    S(0,jj+1) = xab(2) * S(0,jj) + di(jj)*S(0,jj-1)
   enddo
  endif
 
  do jj=1,j
   S(1,jj) = xab(1) * S(0,jj) + di(jj) * S(0,jj-1)
   do ii=2,i
    S(ii,jj) = xab(1) * 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 = dsqrt(pi*inv_p(2))*c*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
! double precision, intent(in)    :: gmma
! 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
 real, intent(in)    :: center(3)
 double precision, intent(in) :: gmma
!DEC$ ATTRIBUTES FORCEINLINE
 integer :: p,q,k
 double precision :: integral
 double precision :: ao_eplf_integral_primitive_oneD
 double precision :: buffer(100)


 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 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 =                              &
   ao_eplf_integral_primitive_oneD(        &
    ao_expo(i,p),                          &
    nucl_coord(ao_nucl(i),1),              &
    ao_power(i,1),                         &
    ao_expo(j,q),                          &
    nucl_coord(ao_nucl(j),1),              &
    ao_power(j,1),                         &
    gmma,                                  &
    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
   ao_eplf_integral += buffer(p)
 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(k,l)
       enddo
     endif
   enddo
 
end function