BEGIN_PROVIDER [ integer, det_num ]
  implicit none
  BEGIN_DOC
! Number of determinants
  END_DOC

 det_num = 1
 call get_determinants_det_num(det_num)
 if (det_num < 1) then
   call abrt(irp_here,'det_num should be > 0')
 endif

END_PROVIDER

BEGIN_PROVIDER  [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
&BEGIN_PROVIDER [ real, det_coef, (det_num) ]
  implicit none

  BEGIN_DOC
! det : Description of the active orbitals of the determinants

! det_coef  : Determinant coefficients
  END_DOC

 if (elec_alpha_num > mo_closed_num) then
  det = 0
  call get_determinants_det_occ(det)
 endif
 det_coef = 0.
 det_coef(1) = 1.
 call get_determinants_det_coef(det_coef)

END_PROVIDER

BEGIN_PROVIDER [ real, mo_occ, (mo_tot_num) ]
 implicit none
 BEGIN_DOC
! Occupation numbers of molecular orbitals
 END_DOC

 call get_mo_basis_mo_occ(mo_occ)

END_PROVIDER



integer function det_exc(k,l,p)
 implicit none
! Degree of excitation+1 between two determinants. Indices are alpha, beta
! The sign is the phase factor

 integer :: k,l,p
 integer :: i, j, jmax

 jmax = elec_num_2(p)-mo_closed_num
 det_exc = 0

 do i=1,jmax
  logical :: found
  found = .False.
  do j=1,jmax
   if (det(j,l,p) == det(i,k,p)) then
     found = .True.
     exit
   endif
  enddo
  if (.not.found) then
    det_exc += 1
  endif

 enddo

 ! Phase

 integer :: nperm
 nperm = 0
 integer, allocatable, save :: buffer(:)
 if (.not. allocated(buffer)) then
   allocate (buffer(0:mo_num-mo_closed_num))
 endif
 do i=1,elec_num_2(p)-mo_closed_num
  buffer(i) = det(i,k,p)
 enddo
 do i=1,elec_num_2(p)-mo_closed_num
  if (buffer(i) /= det(i,l,p)) then
   integer :: m
   m=elec_num_2(p)-mo_closed_num
   do j=i+1,elec_num_2(p)-mo_closed_num
    if (buffer(i) == det(j,l,p)) then  ! found
     m=j
     exit
    endif 
   enddo
   buffer(0) = buffer(i)
   buffer(i) = det(m,l,p)
   buffer(m) = buffer(0)
   nperm += m-i
  endif
 enddo
 det_exc += 1
 det_exc *= (1-2*mod( nperm, 2 ))

end

subroutine get_single_excitation(k,l,m,n,p)
 implicit none
 integer, intent(in)  :: k, l  ! determinant indices
 integer, intent(out) :: m, n  ! m->n excitation
 integer, intent(in)  :: p     ! spin
 
 logical :: found
 integer :: i,j
 m=0
 n=0
 do j=1,elec_num_2(p)-mo_closed_num
  found = .False.
  do i=1,elec_num_2(p)-mo_closed_num
   if (det(j,k,p) == det(i,l,p)) then
    found = .True.
    exit
   endif
  enddo
  if (.not.found) then
    m = det(j,k,p)
    exit
  endif
 enddo

 do i=1,elec_num_2(p)-mo_closed_num
  found = .False.
  do j=1,elec_num_2(p)-mo_closed_num
   if (det(j,k,p) == det(i,l,p)) then
    found = .True.
    exit
   endif
  enddo
  if (.not.found) then
    n = det(i,l,p)
    exit
  endif
 enddo

end

subroutine get_double_excitation(k,l,m,n,r,s,p)
 implicit none
 integer, intent(in)  :: k, l  ! determinant indices
 integer, intent(out) :: m, n  ! m->n excitation
 integer, intent(out) :: r, s  ! r->s excitation
 integer, intent(in)  :: p     ! spin
 
 logical :: found
 integer :: i,j
 m=0
 n=0
 r=0
 s=0
 do j=1,elec_num_2(p)-mo_closed_num
  found = .False.
  do i=1,elec_num_2(p)-mo_closed_num
   if (det(j,k,p) == det(i,l,p)) then
    found = .True.
    exit
   endif
  enddo
  if (.not.found) then
   if (m == 0) then
    m = det(j,k,p)
   else
    r = det(j,k,p)
    exit
   endif
  endif
 enddo

 do i=1,elec_num_2(p)-mo_closed_num
  found = .False.
  do j=1,elec_num_2(p)-mo_closed_num
   if (det(j,k,p) == det(i,l,p)) then
    found = .True.
    exit
   endif
  enddo
  if (.not.found) then
   if (n == 0) then
    n = det(i,l,p)
   else
    s = det(i,l,p)
    exit
   endif
  endif
 enddo

end

BEGIN_PROVIDER [ real, ci_mo, (mo_num,mo_num,3) ]
 implicit none
 BEGIN_DOC  
! Spin Density matrix in the AO basis
 END_DOC
 integer :: i,j,k,l,m,ispin, ik,il
 do ispin=1,3

  do j=1,mo_num
   do i=1,mo_num
    ci_mo(i,j,ispin) = 0.
   enddo
  enddo

  do l=1,det_num
   do m=1,det_num
    real :: factor
    factor = 2.*det_coef(l)*det_coef(m) 
    do il=1,mo_closed_num
     do ik=1,mo_closed_num
      ci_mo(ik,il,ispin) += factor
     enddo
    enddo
   enddo
  enddo

 enddo


 do l=1,det_num
  do m=1,det_num
   factor = det_coef(l)*det_coef(m) 
   do ispin=1,2
    do j=mo_closed_num+1,elec_num_2(ispin)
     ik = det(j-mo_closed_num,l,ispin)
     do il=1,mo_closed_num
      ci_mo(ik,il,ispin) += factor
      ci_mo(il,ik,ispin) += factor
     enddo
     do i=mo_closed_num+1,elec_num_2(ispin)
      il = det(i-mo_closed_num,m,ispin)
      ci_mo(ik,il,ispin) += factor
      ci_mo(il,ik,ispin) += factor
     enddo
    enddo
   enddo

   ispin=3
   do j=mo_closed_num+1,elec_num_2(1)
    ik = det(j-mo_closed_num,l,1)
    do il=1,mo_closed_num
     ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m) 
     ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m) 
    enddo
    do i=mo_closed_num+1,elec_num_2(2)
     il = det(i-mo_closed_num,m,2)
     ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m) 
     ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m) 
    enddo
   enddo

  enddo
 enddo

END_PROVIDER


BEGIN_PROVIDER [ integer, two_e_density_num_max ]
 implicit none
 BEGIN_DOC  
! Number of factors containing the Slater rules 
 END_DOC

 two_e_density_num_max = 2*mo_num

 integer :: k,l
 integer :: exc(3), nact, nact2, p, p2
 integer :: det_exc
 do k=1,det_num
  do l=k,det_num
   exc(1) = abs(det_exc(k,l,1))-1
   exc(2) = abs(det_exc(k,l,2))-1
   exc(3) = exc(1)+exc(2)

   do p=1,2
     p2 = 1+mod(p,2)
     nact  = elec_num_2(p) -mo_closed_num
     nact2 = elec_num_2(p2)-mo_closed_num
     if ( exc(3) == 0 ) then
      two_e_density_num_max += 2*nact*mo_num
     else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
      two_e_density_num_max += 2*mo_num
     else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
      two_e_density_num_max += 2
     else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
      two_e_density_num_max += 1
     endif
   enddo

  enddo
 enddo
 
END_PROVIDER

 BEGIN_PROVIDER [ integer, two_e_density_indice, (4,two_e_density_num_max) ]
&BEGIN_PROVIDER [ real, two_e_density_value, (2,two_e_density_num_max) ]
&BEGIN_PROVIDER [ integer, two_e_density_num ]
 implicit none
 BEGIN_DOC  
! Compact representation of eplf factors
 END_DOC

 integer :: i,j,k,l,m

 integer :: n,p,p2,q
 integer :: ik,il,jk,jl, idx(4)
 real :: phase
 integer :: exc(4), nact, nact2
 real :: det_kl
 integer :: det_exc

 two_e_density_num = 0

 PROVIDE det

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

   exc(1) = det_exc(k,l,1)
   exc(2) = det_exc(k,l,2)
   exc(4) = exc(1)*exc(2)
   exc(1) = abs(exc(1))-1
   exc(2) = abs(exc(2))-1
   exc(3) = exc(1)+exc(2)
   exc(4) = exc(4)/abs(exc(4))
   phase = dble(exc(4))

   det_kl = phase*det_coef(k)*det_coef(l)
   if (k /= l) then
    det_kl += det_kl
   endif

   logical :: notfound
BEGIN_SHELL [ /usr/bin/python ]
code = """
        notfound = .True.
        idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
        do q=1,two_e_density_num
         if (sum(abs(two_e_density_indice(:,q)-idx))) then
          two_e_density_value(1,q) += det_kl
          two_e_density_value(2,q) += det_kl
          notfound = .False.
          exit
         endif
        enddo
        if (notfound) then
         two_e_density_num += 1
         two_e_density_indice(:,two_e_density_num)=idx
         two_e_density_value(1,two_e_density_num) = det_kl
         two_e_density_value(2,two_e_density_num) = det_kl
        endif

        notfound = .True.
        idx = (/ min(%(I)s,%(K)s), max(%(I)s,%(K)s), min(%(J)s,%(L)s), max(%(J)s,%(L)s) /)
        do q=1,two_e_density_num
         if (sum(abs(two_e_density_indice(:,q)-idx))) then
          two_e_density_value(1,q) -= det_kl
          notfound = .False.
          exit
         endif
        enddo
        if (notfound) then
         two_e_density_num += 1
         two_e_density_indice(:,two_e_density_num)=idx
         two_e_density_value(1,two_e_density_num) = -det_kl
         two_e_density_value(2,two_e_density_num) = 0.
        endif
"""

code1 = """
        idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
        notfound = .True.
        do q=1,two_e_density_num
         if (sum(abs(two_e_density_indice(:,q)-idx))) then
          two_e_density_value(1,q) += det_kl
          notfound = .False.
          exit
         endif
        enddo
        if (notfound) then
         two_e_density_num += 1
         two_e_density_indice(:,two_e_density_num)=idx
         two_e_density_value(1,two_e_density_num) = det_kl
         two_e_density_value(2,two_e_density_num) = 0.
        endif

        notfound = .True.
        idx = (/ min(%(I)s,%(K)s), max(%(I)s,%(K)s), min(%(J)s,%(L)s), max(%(J)s,%(L)s) /)
        do q=1,two_e_density_num
         if (sum(abs(two_e_density_indice(:,q)-idx))) then
          two_e_density_value(1,q) -= det_kl
          notfound = .False.
          exit
         endif
        enddo
        if (notfound) then
         two_e_density_num += 1
         two_e_density_indice(:,two_e_density_num)=idx
         two_e_density_value(1,two_e_density_num) = -det_kl
         two_e_density_value(2,two_e_density_num) = 0.
        endif
"""

code2 = """
        notfound = .True.
        idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
        do q=1,two_e_density_num
         if (sum(abs(two_e_density_indice(:,q)-idx))) then
          two_e_density_value(2,q) += det_kl
          notfound = .False.
          exit
         endif
        enddo
        if (notfound) then
         two_e_density_num += 1
         two_e_density_indice(:,two_e_density_num)=idx
         two_e_density_value(1,two_e_density_num) = 0.
         two_e_density_value(2,two_e_density_num) = det_kl
        endif
"""

rep = { \
'CLOSED'          : code%{ 'I':'ik', 'J':'il', 'K':'j', 'L':'j' },
'OPEN_CLOSED'     : code%{ 'I':'j', 'J':'j', 'K':'ik', 'L':'il' },
'OPEN_OPEN_1'     : code1%{ 'I':'ik', 'J':'il', 'K':'jk', 'L':'jl' },
'OPEN_OPEN_2'     : code2%{ 'I':'ik', 'J':'il', 'K':'jk', 'L':'jl' }
}


print """

   do p=1,2
     p2 = 1+mod(p,2)
     nact  = elec_num_2(p) -mo_closed_num
     nact2 = elec_num_2(p2)-mo_closed_num

     if ( exc(3) == 0 ) then
      do n=1,nact
       ik = det(n,k,p)
       il = det(n,l,p)
       do j=1,mo_closed_num
        ! Closed-open shell interactions
        %(CLOSED)s
        !- Open-closed shell interactions
        %(OPEN_CLOSED)s
       enddo

       !- Open-open shell interactions
       do m=1,nact
         jk = det(m,k,p)
         jl = det(m,l,p)
         %(OPEN_OPEN_1)s
       enddo
       do m=1,nact2
         jk = det(m,k,p2)
         jl = det(m,l,p2)
         %(OPEN_OPEN_2)s
       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)

      do j=1,mo_closed_num
       !- Open-closed shell interactions
       %(CLOSED)s
       !- Closed-open shell interactions
       %(OPEN_CLOSED)s
      enddo

      !- Open-open shell interactions
      do m=1,nact
       jk = det(m,k,p)
       jl = det(m,l,p)
       %(OPEN_OPEN_1)s
      enddo
      do m=1,nact2
       jk = det(m,k,p2)
       jl = det(m,l,p2)
       %(OPEN_OPEN_2)s
      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)
      %(OPEN_OPEN_1)s

     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)
      %(OPEN_OPEN_2)s
    
     endif

   enddo
"""%(rep)
END_SHELL

  enddo
 enddo
 
END_PROVIDER

BEGIN_PROVIDER [ real, one_e_density_mo, (mo_active_num,mo_active_num,2) ]
 implicit none
 BEGIN_DOC  
! One electron spin density matrix in MO space
 END_DOC
 integer :: i,j,k,l,p, il, jl
 do p=1,2
  do i=1,mo_active_num
   do j=1,mo_active_num
    one_e_density_mo(j,i,p) = 0.
   enddo
  enddo
 enddo
 real :: ckl, phase
 integer :: exc(4), det_exc
 do k=1,det_num
  do l=k,det_num
   exc(1) = det_exc(k,l,1)
   exc(2) = det_exc(k,l,2)
   exc(4) = exc(1)*exc(2)
   exc(1) = abs(exc(1))-1
   exc(2) = abs(exc(2))-1
   exc(3) = exc(1)+exc(2)
   exc(4) = exc(4)/abs(exc(4))
   phase = dble(exc(4))
   ckl = det_coef(k)*det_coef(l)*phase
   do p=1,2
    if (exc(3) == 0) then
      do i=1,elec_num_2(p)-mo_closed_num
       il = det(i,k,p) - mo_closed_num
       one_e_density_mo(il,il,p) += ckl
      enddo
    else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
      call get_single_excitation(k,l,il,jl,p)
      jl -= mo_closed_num
      il -= mo_closed_num
      one_e_density_mo(il,jl,p) += ckl
      one_e_density_mo(jl,il,p) += ckl
    endif
   enddo
  enddo
 enddo

END_PROVIDER