subroutine order_quad_indices_2_2(ij_ab, kl_cd, phase, ijab_klcd)
integer, intent(inout) :: ij_ab(2, 2), kl_cd(2, 2)
real(dp), intent(out) :: phase
integer, intent(out) :: ijab_klcd(8)
character(*), parameter :: this_routine = "order_quad_indices_2_2"
integer :: ij(2), ab(2), kl(2), cd(2)
integer :: n
! the correctly order t_ij^ab * t_kl^cd has the sign convention -1!
! and we can be sure that the inputed (ij),(ab) etc. are ordered i < j
! within the T2 operators, since only those get stored by convention
! and also that i < k and a < c
ij = ij_ab(1, :)
ab = ij_ab(2, :)
kl = kl_cd(1, :)
cd = kl_cd(2, :)
n = 1
! assert i < k and a < c
ASSERT(ij_ab(1, 1) < kl_cd(1, 1))
ASSERT(ij(1) < ij(2))
ASSERT(cd(1) < cd(2))
ASSERT(ab(1) < ab(2))
ASSERT(cd(1) < cd(2))
! with this: i am not so sure: a < c ..
! ok i can not be sure that a < c! so i have to sort more here!
! sort the electrons:
! if j < k everything is fine and in order!
if (ij_ab(1, 2) < kl_cd(1, 1)) then
! everything is fine..
else if (ij_ab(1, 2) < kl_cd(1, 2)) then
! if j < l then i have to switch j and k
call swap(ij(2), kl(1))
n = n + 1
else
! this means j > l so we have to swap j -> l and then k -> l
call swap(ij(2), kl(2))
call swap(ij(2), kl(1))
! wich leaves the phase unchanged
end if
! and check if we did everything correctly:
ASSERT(ij(1) < ij(2))
ASSERT(kl(1) < kl(2))
ASSERT(ij(2) < kl(1))
! i am not so sure any more if a < c anymore
! now sort the orbitals but here a > c is also possible!
if (ij_ab(2, 1) < kl_cd(2, 1)) then
! so a < c
if (ij_ab(2, 2) < kl_cd(2, 1)) then
! b < c so in this case nothing has to be done!
else if (ij_ab(2, 2) < kl_cd(2, 2)) then
! b < d: so b and c have to be swapped!
call swap(ab(2), cd(1))
n = n + 1
else
! b > d: so switch b -> d and d -> c
call swap(ab(2), cd(2))
call swap(ab(2), cd(1))
! and this does not change the phase
end if
else
! so this means a > c which implies b > c also!
! so if b < d we know everything
if (ij_ab(2, 2) < kl_cd(2, 2)) then
! then c < a < b < d
call swap(ab(1), cd(1))
call swap(ab(2), cd(1))
! so no phase factor!
else if (ij_ab(2, 1) > kl_cd(2, 2)) then
! here b and a > d
! so c < d < a < b
call swap(ab(1), cd(1))
call swap(ab(2), cd(2))
! also here no phase factor!
else
! this means c < a < d < b
call swap(ab(1), ab(2))
! ba, cd
call swap(cd(1), cd(2))
! ba, dc
call swap(ab(1), cd(2))
! cadb
! here we get a phase
n = n + 1
end if
end if
ASSERT(ab(1) < ab(2))
ASSERT(cd(1) < cd(2))
ASSERT(ab(2) < cd(1))
! now switch the indices
ij_ab(1, :) = ij
ij_ab(2, :) = ab
kl_cd(1, :) = kl
kl_cd(2, :) = cd
! and also store a linear index
ijab_klcd = [ij, ab, kl, cd]
! and calculate the appropriate phase
phase = (-1.0_dp)**n
end subroutine order_quad_indices_2_2