4 symmetries

a.f90

integer recursive function ijkl2intindex2(i, j, k, l) result(r)
! Index into the int2(:) array using the chemistry notation (ij|kl)
! Only the 4 general symmetries are assumed.
! Here is how to iterate over all (ij|kl) combinations:
!    ijkl = 1
!    do i = 1, n
!        do j = 1, i
!            do k = 1, i
!                do l = 1, i
!                    if (i == j .and. k < l) cycle
!                    if (i == k .and. j < l) cycle
!                    if (i == l .and. j < k) cycle
!                    call assert(ijkl2intindex2(i, j, k, l) == ijkl)
!                    ijkl = ijkl + 1
!                end do
!            end do
!        end do
!    end do
integer, intent(in) :: i, j, k, l
if (i < j .or. (i == j .and. k < l)) then
    r = ijkl2intindex2(j, i, l, k)
else if (i < k .or. (i == k .and. j < l)) then
    r = ijkl2intindex2(k, l, i, j)
else if (i < l .or. (i == l .and. j < k)) then
    r = ijkl2intindex2(l, k, j, i)
else
    r = getni(i) + getnj(i, j) + getnk(i, j, k) + getnl(i, j, k, l)
end if
end function

integer pure function getnl(i, j, k, l) result(n)
integer, intent(in) :: i, j, k, l
n = l
if (i == j) n = min(k, n)
if (i == k) n = min(j, n)
if (i <= n .and. j < k) n = n - 1
end function

integer pure function getnk(i, j, k_) result(n)
integer, intent(in) :: i, j, k_
integer :: k
n = 0
do k = 1, k_ - 1
    n = n + getnl(i, j, k, i)
end do
end function

integer pure function getnj(i, j_) result(n)
integer, intent(in) :: i, j_
integer :: j
n = 0
do j = 1, j_ - 1
    n = n + getnk(i, j, i+1)
end do
end function

integer pure function getni(i) result(n)
integer, intent(in) :: i
n = (i-1)**2 * ((i-1)**2 + 3) / 4
end function