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

 integer :: dl(mo_closed_num), dk(mo_closed_num), buffer(0:mo_closed_num)
 do i=1,jmax
  dk(i) = det(i,k,p)
  dl(i) = det(i,l,p)
  buffer(i) = dk(i)
 enddo

 integer :: kmax
 logical :: notfound
 do i=1,jmax
  notfound = .True.
  do j=1,jmax
   notfound = notfound .and. (dl(j) /= dk(i))
  enddo
  if (notfound) then
    det_exc += 1
  endif
 enddo

 ! Phase

 integer :: nperm
 nperm = 0
 do i=1,jmax
  if (buffer(i) /= dl(i)) then
   integer :: m
   m=jmax
   do j=i+1,jmax
    if (buffer(i) == dl(j)) then  ! found
     m=j
     exit
    endif 
   enddo
   buffer(0) = buffer(i)
   buffer(i) = dl(m)
   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 :: notfound
 integer :: i,j
 m=0
 n=0

 integer :: dl(elec_alpha_num-mo_closed_num), dk(elec_alpha_num-mo_closed_num), buffer(0:elec_alpha_num-mo_closed_num)
 integer :: jmax
 jmax = elec_num_2(p)-mo_closed_num

 do i=1,jmax
  dk(i) = det(i,k,p)
  dl(i) = det(i,l,p)
 enddo

 do j=1,jmax
  notfound = .True.
  do i=1,jmax
   notfound = notfound .and. (dk(j) /= dl(i))
  enddo
  if (notfound) then
    m = dk(j)
    exit
  endif
 enddo

 do i=1,jmax
  notfound = .True.
  do j=1,jmax
   notfound = notfound .and. (dk(j) /= dl(i))
  enddo
  if (notfound) 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 :: notfound
 integer :: i,j
 m=0
 n=0
 r=0
 s=0

 integer :: dl(elec_alpha_num-mo_closed_num), dk(elec_alpha_num-mo_closed_num), buffer(0:elec_alpha_num-mo_closed_num)
 integer :: jmax
 jmax = elec_num_2(p)-mo_closed_num

 do i=1,jmax
  dk(i) = det(i,k,p)
  dl(i) = det(i,l,p)
 enddo

 do j=1,jmax
  notfound = .True.
  do i=1,jmax
   notfound = notfound .and. (dk(j) /= dl(i))
  enddo
  if (notfound) then
   if (m == 0) then
    m = dk(j)
   else
    r = dk(j)
    exit
   endif
  endif
 enddo

 do j=1,jmax
  notfound = .True.
  do i=1,jmax
   notfound = notfound .and. (dk(j) /= dl(i))
  enddo
  if (notfound) then
   if (n == 0) then
    n = dl(j)
   else
    s = dl(j)
    exit
   endif
  endif
 enddo

end


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 = 0
 call get_density_matrix_two_num(two_e_density_num_max)
 if (two_e_density_num_max /= 0) then
   return
 endif

 two_e_density_num_max = 2*mo_num
 integer :: k,l
 integer :: exc(3), nact, p
 integer :: det_exc
 do k=1,det_num
  if (abs(det_coef(k)) < 1.e-5) then
    cycle
  endif
  do l=k,det_num
   if ( (k /= l).and.(abs(det_coef(k)*det_coef(l)) < 1.e-5) ) then
     cycle
   endif
   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
     nact  = elec_num_2(p) -mo_closed_num
     if ( exc(3) == 0 ) then
      two_e_density_num_max += nact*(mo_closed_num*4+nact*3)
     else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
      two_e_density_num_max += mo_closed_num*4+nact*3
     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
 
 call set_density_matrix_two_num(two_e_density_num_max)
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

 two_e_density_indice(1,1) = -1
 call get_density_matrix_two_indice(two_e_density_indice)
 call get_density_matrix_two_value(two_e_density_value)
 call get_density_matrix_two_num(two_e_density_num)
 if (two_e_density_indice(1,1) /= -1) then
   return
 endif

 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, det_k
 integer :: det_exc

 two_e_density_num = 0

 PROVIDE det elec_alpha_num

 do k=1,det_num
  det_k = det_coef(k)
  if ( abs(det_k) < 1.e-5) then
    cycle
  endif
  do l=k,det_num
   det_kl = det_k*det_coef(l)
   if ( (k /= l).or.(abs(det_kl) < 1.e-5) ) then
     cycle
   endif

   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
   if (k /= l) then
    det_kl += det_kl
   endif

   logical :: notfound
BEGIN_SHELL [ /usr/bin/python ]
code = """
        notfound = .True.
        idx(1) = min(%(I)s,%(J)s)
        idx(2) = max(%(I)s,%(J)s)
        idx(3) = min(%(K)s,%(L)s)
        idx(4) = max(%(K)s,%(L)s)
        do q=1,two_e_density_num
         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
              (two_e_density_indice(2,q)==idx(2)) .and. &
              (two_e_density_indice(3,q)==idx(3)) .and. &
              (two_e_density_indice(4,q)==idx(4)) ) 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(1,two_e_density_num)=idx(1)
         two_e_density_indice(2,two_e_density_num)=idx(2)
         two_e_density_indice(3,two_e_density_num)=idx(3)
         two_e_density_indice(4,two_e_density_num)=idx(4)
         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(1) = min(%(I)s,%(K)s)
        idx(2) = max(%(I)s,%(K)s)
        idx(3) = min(%(J)s,%(L)s)
        idx(4) = max(%(J)s,%(L)s)
        do q=1,two_e_density_num
         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
              (two_e_density_indice(2,q)==idx(2)) .and. &
              (two_e_density_indice(3,q)==idx(3)) .and. &
              (two_e_density_indice(4,q)==idx(4)) ) 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(1,two_e_density_num)=idx(1)
         two_e_density_indice(2,two_e_density_num)=idx(2)
         two_e_density_indice(3,two_e_density_num)=idx(3)
         two_e_density_indice(4,two_e_density_num)=idx(4)
         two_e_density_value(1,two_e_density_num) = -det_kl
         two_e_density_value(2,two_e_density_num) = 0.
        endif
"""

code1 = """
        idx(1) = min(%(I)s,%(J)s)
        idx(2) = max(%(I)s,%(J)s)
        idx(3) = min(%(K)s,%(L)s)
        idx(4) = max(%(K)s,%(L)s)
        notfound = .True.
        do q=1,two_e_density_num
         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
              (two_e_density_indice(2,q)==idx(2)) .and. &
              (two_e_density_indice(3,q)==idx(3)) .and. &
              (two_e_density_indice(4,q)==idx(4)) ) 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(1,two_e_density_num)=idx(1)
         two_e_density_indice(2,two_e_density_num)=idx(2)
         two_e_density_indice(3,two_e_density_num)=idx(3)
         two_e_density_indice(4,two_e_density_num)=idx(4)
         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(1) = min(%(I)s,%(K)s)
        idx(2) = max(%(I)s,%(K)s)
        idx(3) = min(%(J)s,%(L)s)
        idx(4) = max(%(J)s,%(L)s)
        do q=1,two_e_density_num
         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
              (two_e_density_indice(2,q)==idx(2)) .and. &
              (two_e_density_indice(3,q)==idx(3)) .and. &
              (two_e_density_indice(4,q)==idx(4)) ) 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(1,two_e_density_num)=idx(1)
         two_e_density_indice(2,two_e_density_num)=idx(2)
         two_e_density_indice(3,two_e_density_num)=idx(3)
         two_e_density_indice(4,two_e_density_num)=idx(4)
         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(1) = min(%(I)s,%(J)s)
        idx(2) = max(%(I)s,%(J)s)
        idx(3) = min(%(K)s,%(L)s)
        idx(4) = max(%(K)s,%(L)s)
        do q=1,two_e_density_num
         if ( (two_e_density_indice(1,q)==idx(1)) .and. &
              (two_e_density_indice(2,q)==idx(2)) .and. &
              (two_e_density_indice(3,q)==idx(3)) .and. &
              (two_e_density_indice(4,q)==idx(4)) ) 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(1,two_e_density_num)=idx(1)
         two_e_density_indice(2,two_e_density_num)=idx(2)
         two_e_density_indice(3,two_e_density_num)=idx(3)
         two_e_density_indice(4,two_e_density_num)=idx(4)
         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
 
 call set_density_matrix_two_num(two_e_density_num)
 call set_density_matrix_two_indice(two_e_density_indice)
 call set_density_matrix_two_value(two_e_density_value)

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

 if (mo_active_num == 0) then
   return
 endif
 one_e_density_mo(1,1,1) = -1.
 call get_density_matrix_one(one_e_density_mo)
 if (one_e_density_mo(1,1,1) /= -1.) then
   return
 endif

 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
 PROVIDE mo_closed_num elec_num_2 elec_alpha_num 
 do k=1,det_num
  do l=k,det_num
   ckl = det_coef(k)*det_coef(l)
   if ( (k /= l).and.(abs(ckl) < 1.e-5) ) then
     cycle
   endif
   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 *= 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

 call set_density_matrix_one(one_e_density_mo)

END_PROVIDER