# How to implement efficient indexing function for two particle integrals <ij|kl>?

This is a simple symmetry enumeration problem. I give the full background here, but no knowledge of quantum chemistry is needed.

The two particle integral $\langle ij|kl\rangle$ is: $$\langle ij|kl\rangle = \int {\psi_i^*({\bf x})\psi_j^*({\bf x}')\psi_k({\bf x})\psi_l({\bf x}') \over | {\bf x} - {\bf x}' |}\, d^3 x \,d^3 x'$$ And it has the following 4 symmetries: $$\langle ij|kl\rangle = \langle ji|lk\rangle = \langle kl|ij\rangle = \langle lk|ji\rangle$$ I have a function that calculates the integral and stores them in a 1D array int2, indexed as follows:

int2(ijkl2intindex2(i, j, k, l))


where the function ijkl2intindex2 returns a unique index, taking the above symmetries into account. The only requirement is that if you loop over all combinations of i, j, k, l (from 1 to n each), it will fill the int2 array consecutively and it will assign the same index to all ijkl combinations that are related by the above 4 symmetries.

My current implementation in Fortran is here. It is very slow. Does anybody know how to do this effectively? (In any language.)

Hint: if the orbitals $\psi_i({\bf x})$ are real, then in addition to the above symmetries, one can exchange $i\leftrightarrow k$ and $j\leftrightarrow l$ so we get 8 symmetries total: $$\langle ij|kl\rangle = \langle ji|lk\rangle = \langle kj|il\rangle = \langle il|kj\rangle =$$ $$= \langle kl|ij\rangle = \langle lk|ji\rangle = \langle il|kj\rangle = \langle kj|il\rangle$$ and then one can implement a very fast function for indexing it, see my implementation here. I would like to find some efficient indexing scheme for the cases when the orbitals are not real.

Note: the functions that I implemented actually accept the four numbers $i$, $j$, $k$, $l$ in a so called "chemistry" notation $(ij|kl) = \langle ik|jl\rangle$, i.e. the $j$ and $k$ arguments are interchanged, but this is not important.

• Off-topic, but am I the only person appalled by the notation $d^3x$? Kill it, kill it with fire! – n00b May 23 '12 at 11:31
• $d^3 x$ is just a commonly used shortcut for $d x_1 d x_2 d x_3$ where ${\bf x} = (x_1, x_2, x_3)$. Instead of writing it all out explicitly, it's much simpler to just use $d^3 x$ and all is clear. – Ondřej Čertík May 23 '12 at 19:03
• I can see clearly what it is: misleading, hideous and in need of destruction. There is no $x$ in $\mathbf{x} =(x_1,x_2,x_3)$ so why is it the integration variable? Why not just $\mathrm{d}\mathbf{x}$ ? I feel ill just looking at it, and now I need a lie down. – n00b May 24 '12 at 8:56
• @n00b, I think that $d^3 x$ is preferred because it also specifies the dimension of the integral (very important, as the integral gives different results in 1D, 2D and 3D). – Ondřej Čertík May 24 '12 at 22:32

[Edit: 4th time's the charm, at last something sensible]

I got to this backwards: I started with another answer showing the filter-based approach, and used that to generate all the valid combinations for a series of values of $n$, and looked the sequence up on the online integer sequence database. The number of combinations is $n^2(n^2+3)$, which looks unlikely (why the 3?). It is also $t(t(n))+t(t(n-1))$ where $t(a)$ is the triangular number of $a$, $t(a)=a(a+1)/2$. Having got this, we need to know why.

The first term is simpler - the pairs of pairs where the $i \ge j$ and $tid(i,j) \ge tid(k,l)$, where $tid(a,b)$ is the triangular index of $a,b$. That is fulfilled by a function like this:

def ascendings(n):
idx = 0
for i in range(1,n+1):
for j in range(1,i+1):
for k in range(1,i):
for l in range(1,k+1):
idx = idx + 1
print(i,j,k,l)
k=i
for l in range(1,j+1):
idx = idx + 1
print(i,j,k,l)
return idx


where the second $l$ loop is because we cannot nest the $l$ loop entirely inside the $k$ loop without an if/skip to check the triangular indices.

The second term $t(t(n-1))$ is where the first pair is ascending, and the second pair is descending (note, no ==, as they are handled above).

def mixcendings(n):
idx = 0
for j in range(2,n+1):
for i in range(1,j):
for k in range(1,j):
for l in range(1,k):
print(i,j,k,l)
idx = idx + 1
k=j
for l in range(1,i+1):
print(i,j,k,l)
idx = idx + 1
return idx


The combination of both of these gives the complete set, so putting both loops together gives us the complete set of indices.

A significant problem is that these patterns are difficult to calculate for an arbitrary i,j,k,l. So I would suggest a map that yields the index given i,j,k,l. Frankly, if you are doing this at all, you might as well use the generate+filter approach, because you only need to do that once for a given $n$. The plus side of the above method is that you at least have a predictable loop structure.

In python we can write the following iterator to give us the idx and i,j,k,l values for each different scenario:

def iterate_quad(n):
idx = 0
for i in range(1,n+1):
for j in range(1,i+1):
for k in range(1,i):
for l in range(1,k+1):
idx = idx + 1
yield (idx,i,j,k,l)
#print(i,j,k,l)
k=i
for l in range(1,j+1):
idx = idx + 1
yield (idx,i,j,k,l)

for i in range(2,n+1):
for j in range(1,i):
for k in range(1,i):
for l in range(1,k):
idx = idx + 1
yield (idx,i,j,k,l)
k=i
for l in range(1,j+1):
idx = idx + 1
yield (idx,i,j,k,l)


In fortran we would just have to run the loop and store the values. We can use a simple index to store the i,j,k,l combination as a single value ($i n^3 + j n^2 + k n + l$), and store these values in an array whose index is the same as the index above. We can then iterate over this array and retrieve the i,j,k,l from the values. To get the idx for arbitrary i,j,k,l would require a reverse map and a filter to handle the symmetry, although we could probably construct a function from the above structure. The idx array generation function in fortran would be:

integer function squareindex(i,j,k,l,n)
integer,intent(in)::i,j,k,l,n
squareindex = (((i-1)*n + (j-1))*n + (k-1))*n + l
end function

integer function generate_order_array(n,arr)
integer,intent(in)::n,arr(*)
integer::total,idx,i,j,k,l
total = n**2 * (n**2 + 3)
reshape(arr,total)
idx = 0
do i=1,n
do j=1,i
do k=1,i-1
do l=1,k
idx = idx+1
arr(idx) = squareindex(i,j,k,l,n)
end do
end do
k=i
do l=1,j
idx = idx+1
arr(idx) = squareindex(i,j,k,l,n)
end do
end do
end do

do i=2,n
do j=1,i-1
do k=1,i-1
do l=1,j
idx = idx+1
arr(idx) = squareindex(i,j,k,l,n)
end do
end do
k=i
do l=1,j
idx = idx+1
arr(idx) = squareindex(i,j,k,l,n)
end do
end do
end do

generate_order_array = idx
end function


And then loop over it thus:

maxidx = generate_order_array(n,arr)
do idx=1,maxidx
i = idx/(n**3) + 1
t_idx = idx - (i-1)*n**3
j = t_idx/(n**2) + 1
t_idx = t_idx - (j-1)*n**2
k = t_idx/n + 1
t_idx = t_idx - (k-1)*n
l = t_idx

! now have i,j,k,l, so do stuff
! ...
end do

• Hi Phil, thanks a lot for the answer! I tested it, and there are two problems. For example idx_all(1, 2, 3, 4, 4) == idx_all(1, 2, 4, 3, 4) = 76. But <12|34> /= <12|43>. It is equal only if the orbitals are real. So your solution seems to be for the case of 8 symmetries (see my Fortran example above for a simpler version, the ijkl2intindex()). The second problem is that the indices are not consecutive, I pasted the results here: gist.github.com/2703756. Here are the correct results from my ijkl2intindex2() routine above: gist.github.com/2703767. – Ondřej Čertík May 15 '12 at 18:03
• @OndřejČertík: You want a sign associated? have idxpair return a sign if you switched order. – Deathbreath May 16 '12 at 15:31
• OndřejČertík: I see the difference now. As @Deathbreath points out, you can negate the index, but that won't be as clean for the overall loop. I'll have a think and update it. – Phil H May 16 '12 at 16:02
• Actually, negating the index won't fully work as the idxpair will get the value wrong. – Phil H May 16 '12 at 16:05
• @PhilH: Don't negate the index $$<ij|kl>=-<ji|kl>=-<ij|lk> =<ji|lk>$$. Just add a return variable sign to idxpair. Then the answer is $$ijkl[idxpair(index_{ij},index_{kl},,)]sign_{ij}sign_{kl}$$. – Deathbreath May 16 '12 at 17:15

Here is an idea of using a simple space filling curve modified to return same key for the symmetry cases (all code snippets are in python).

# Simple space-filling curve
def forge_key(i, j, k, l, n):
return i + j*n + k*n**2 + l*n**3

# Considers the possible symmetries of a key
def forge_key_symmetry(i, j, k, l, n):
return min(forge_key(i, j, k, l, n),
forge_key(j, i, l, k, n),
forge_key(k, l, i, j, n),
forge_key(l, k, j, i, n))


Notes:

• The example is python but if you inline the functions into your fortran code and unroll your inner loop for (i, j, k, l), you should get decent performance.
• You could compute the key using floats and then convert the key into integer to use as an index, this would allow the compiler to use the floating point units (e.g. AVX is available).
• If N is a power of 2 then the multiplications would be just bit shifts.
• The treatment for the symmetries is not efficient in memory (i.e. it does not produce a continuous indexing) and uses around of 1/4 of the total index array entries.

Here is a test example for n = 2:

for i in range(n):
for j in range(n):
for k in range(n):
for l in range(n):
key = forge_key_symmetry(i, j, k, l, n)
print i, j, k , l, key


Output for n = 2:

i j k l key
0 0 0 0 0
0 0 0 1 1
0 0 1 0 1
0 0 1 1 3
0 1 0 0 1
0 1 0 1 5
0 1 1 0 6
0 1 1 1 7
1 0 0 0 1
1 0 0 1 6
1 0 1 0 5
1 0 1 1 7
1 1 0 0 3
1 1 0 1 7
1 1 1 0 7
1 1 1 1 15


If of interest, the inverse function of forge_key is:

# Inverse of forge_key
def split_key(key, n):
d = key / n**3
c = (key - d*n**3) / n**2
b = (key - c*n**2 - d*n**3) / n
a = (key - b*n - c*n**2 - d*n**3)
return (a, b, c, d)

• Did you mean "if n is a power of 2" instead of multiple of 2? – Aron Ahmadia May 17 '12 at 12:22
• yes, thanks Aron. I wrote this answer just before going for dinner and Hulk was writing. – fcruz May 17 '12 at 14:11
• Clever! However, isn't the maximum index n^4 (or n^4-1 if you start from 0)? The problem is that for the basis size that I want to be able to do, it won't fit into memory. With consecutive index, the size of the array is n^2 * (n^2 + 3) / 4. Hm, that's only about 1/4 of the full size anyway. So maybe I should not worry about the factor of 4 in memory consumption. Still though, there must be some way to encode the correct consecutive index using only these 4 symmetries (better than my ugly solution in my post, where I need to do double loops). – Ondřej Čertík May 17 '12 at 22:56
• yeah, that is right! I don't know how to elegantly solve (without sorting and renumbering) the index but the leading term in memory usage is O(N^4). The factor of 4 should make a small difference in memory for large N. – fcruz May 18 '12 at 9:05

Is this not just the generalization of the packed symmetric matrix indexing problem? The solution there is offset(i,j) = i*(i+1)/2+j, is it not? Can't you double down on this and index a doubly-symmetric 4D array? The implementation requiring branching seems unnecessary.