BEGIN_PROVIDER  [ real, eplf_gamma ]
 implicit none
 BEGIN_DOC  
! Value of the gaussian for the EPLF
 END_DOC
 eplf_gamma = 1000.
END_PROVIDER

BEGIN_PROVIDER [ double precision, 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 :: eplf_integral
 do i=1,ao_num
  do j=i,ao_num
    eplf_integral_matrix(j,i) = eplf_integral(j,i,eplf_gamma,point)
    eplf_integral_matrix(i,j) = 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, k, l, m, n
 double precision :: thr

 thr = 1.d-12 / eplf_gamma
 eplf_up_up = 0.
 eplf_up_dn = 0.

 PROVIDE mo_coef_transp

 do m=1,ao_num
  do n=1,ao_num
   double precision :: ao_mn
   ao_mn = eplf_integral_matrix(m,n)

   if (abs(ao_mn) > thr) then
    
    do k=1,ao_num
     do l=1,ao_num
      double precision :: ao_klmn
      ao_klmn = ao_mn*ao_value_p(k)*ao_value_p(l)
    
      if (abs(ao_klmn) > thr) then

       do j=1,elec_alpha_num
        if (mo_coef_transp(j,n) /= 0.d0) then
          do i=1,elec_alpha_num
           eplf_up_up = eplf_up_up + ao_klmn* &
             mo_coef_transp(i,k)*mo_coef_transp(j,n) * &
            (mo_coef_transp(i,l)*mo_coef_transp(j,m) - &
             mo_coef_transp(i,m)*mo_coef_transp(j,l) )
          enddo
        endif
       enddo

       do j=1,elec_beta_num
        if (mo_coef_transp(j,n) /= 0.d0) then
          do i=1,elec_beta_num
           eplf_up_up = eplf_up_up + ao_klmn* &
             mo_coef_transp(i,k)*mo_coef_transp(j,n) * &
            (mo_coef_transp(i,l)*mo_coef_transp(j,m) - &
             mo_coef_transp(i,m)*mo_coef_transp(j,l) )
          enddo
        endif
       enddo
    
       do j=1,elec_beta_num
        if ( (mo_coef_transp(j,n) /= 0.d0).and. &
             (mo_coef_transp(j,m) /= 0.d0) ) then
          do i=1,elec_alpha_num
           eplf_up_dn = eplf_up_dn + 2.d0*ao_klmn* &
            mo_coef_transp(i,k)*mo_coef_transp(j,n) * &
            mo_coef_transp(i,l)*mo_coef_transp(j,m)
          enddo
        endif
       enddo

      endif
    
     enddo
    enddo

   endif

  enddo
 enddo

END_PROVIDER


BEGIN_PROVIDER [ real, eplf_value ]
 implicit none
 BEGIN_DOC  
! Value of the EPLF at the current point.
 END_DOC
 double precision :: aa, ab

 aa = eplf_up_up
 ab = eplf_up_dn
 aa = min(1.d0,aa)
 ab = min(1.d0,ab)
 aa = max(1.d-30,aa)
 ab = max(1.d-30,ab)
 aa = -dlog(aa)/eplf_gamma
 ab = -dlog(ab)/eplf_gamma
 aa = dsqrt(aa)
 ab = dsqrt(ab)

 eplf_value = (aa-ab)/(aa+ab)

END_PROVIDER


double precision function 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=1000
 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 = dtemp + &
   (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2))
  x = x+dx
 enddo
 eplf_integral_primitive_oneD_numeric = dtemp*dx

end function

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

  do q=1,ao_prim_num(j)
   do p=1,ao_prim_num(i)
    integral(p,q) =                         &
    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))   *                         &
    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))   *                         &
    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))
   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
 
  eplf_integral_numeric = 0.
  do q=1,ao_prim_num(j)
   do p=1,ao_prim_num(i)
    eplf_integral_numeric = eplf_integral_numeric + integral(p,q)
   enddo
  enddo
  
end function
 
double precision function 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 :: S(0:i+1,0:j)
 double precision :: inv_p, di(max(i,j)), dj(j), c, c1

 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
 S(0,0) = dsqrt(pi*inv_p)*c*c1
 
 ! 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

 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

 eplf_integral_primitive_oneD = S(i,j)

end function

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

 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) =                         &
   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))   *                         &
   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))   *                         &
   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))
  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

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