8 symmetries

a.f90

integer pure function ijkl2intindex(i_, j_, k_, l_) result(r)
! Index into the int2(:) array using the chemistry notation (ij|kl)
! Here is how to iterate over all (ij|kl) combinations:
!    ijkl = 1
!    do i = 1, n
!        do j = 1, i
!            do k = 1, n
!                do l = 1, k
!                    if ((i-1)*i/2 + j < (k-1)*k/2 + l) cycle
!                    call assert(ijkl2intindex(i, j, k, l) == ijkl)
!                    ijkl = ijkl + 1
!                end do
!            end do
!        end do
!    end do
integer, intent(in) :: i_, j_, k_, l_
integer :: i, j, k, l
integer :: ij, kl
if (i_ < j_) then
    i = j_
    j = i_
else
    i = i_
    j = j_
end if
if (k_ < l_) then
    k = l_
    l = k_
else
    k = k_
    l = l_
end if
ij = (i-1)*i/2 + j
kl = (k-1)*k/2 + l
if (ij < kl) then
    r = (kl-1)*kl/2 + ij
else
    r = (ij-1)*ij/2 + kl
end if
end function