Skip to content

Instantly share code, notes, and snippets.

@certik

certik/a.f90

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