pure function oldLMatInd(aI, bI, cI, iI, jI, kI) result(index)
! Indexing function with a 12-fold symmetry: symmetric with respect to
! permuting (a,i), (b,j) and (c,k) and with exchange (a,b,c)<->(i,j,k)
! Input: a,b,c - orbital indices of electrons
! i,j,k - orbital indices of holes
! Output: index - contiguous index I(a,b,c,i,j,k) with the aforementioned symmetry
integer(int64), value, intent(in) :: aI, bI, cI ! occupied orb indices
integer(int64), value, intent(in) :: iI, jI, kI ! unoccupied orb
integer(int64) :: index
integer(int64) :: a, b, c, i, j, k
! Somehow, I cannot directly use aI-kI, even though they are passed by value
a = aI
b = bI
c = cI
i = iI
j = jI
k = kI
! we store the permutation where a < b < c (regardless of i,j,k)
! or i < j < k, depending on (permuted) a < i
! sort such that the ordered indices start with the smallest index
if (minval((/a, b, c/)) > minval((/i, j, k/))) then
call swap(a, i)
call swap(b, j)
call swap(c, k)
end if
! -> create the ordered permutation on ap,bp,cp
call sort2Els(a, b, i, j)
call sort2Els(b, c, j, k)
call sort2Els(a, b, i, j)
! indexing function: there are three ordered indices (ap,bp,cp)
! and three larger indices (ip,jp,kp)
! the last larger index kp, it is the contigous index, then follow (jp,cp)
! then (ip,bp) and then the smallest index ap
index = k + nBI * (j - 1) + nBI**2 * (i - 1) + nBI**3 * (a - 1) + nbI**3 * (b - 1) * b / 2 + &
nBI**3 * (c + 1) * (c - 1) * c / 6
contains
! sorts the indices a,b and i,j with respect to the
! ordering selected in iPermute
pure subroutine sort2Els(r, s, p, q)
integer(int64), intent(inout) :: r, s, p, q
if (r > s) then
call swap(r, s)
call swap(p, q)
end if
end subroutine sort2Els
end function oldLMatInd