Skip to content

Instantly share code, notes, and snippets.

@certik

certik/a.f90

Created May 14, 2012
Embed
What would you like to do?
4 symmetries
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment