Two electron integrals are of the form (ij|kl), with permutation rules for the 4 indices. It is not possible to store all the (ij|kl) integrals in memory under the form of a 4-dimensional array, so we store instead all the non-zero integrals in a data structure, similar to a hash table.

We do not use hash tables because :

Rehashing the table as it grows is too costly
We need some locality if possible : (ij|kl) should be close to (ij|k l+1)
We tried boost and C++/STL unordered maps, and the hash tables of Google but it was still too slow

We do not use binary search trees because:

We need to be able to add/update some data in a shared-memory parallel environment
It makes lots of jumps among memory pages before the data is found

So we decided to implement our own key/value data structure well adapted to our problem. This is implemented in the Utils/map_module.f90 file.

I would like to re-implement this data-structure in C as an open-source independent library, so that it could be easily used in many other programs and other languages than Fortran. Also, with a library one has the freedom to add some intrinsics for a better use of the hardware. If you are willing to contribute, here is the address of this project : https://github.com/scemama/sparse_array

The Hash function

The key is a function that maps the four indices i,j,k,l into a unique 64-bit integer, and which returns the same value when then permutation symmetry operations are applied to the indices. As it is a tiny function, it is written in the module that uses it intensively such that it can be inlined by the compiler, namely the Integrals_Bielec/mo_bi_integrals.irp.f file.

subroutine mo_bielec_integrals_index(i,j,k,l,i1)
  implicit none
  integer, intent(in) :: i,j,k,l
  integer(key_kind), intent(out) :: i1
  integer(key_kind) :: p,q,r,s,i2
  p = min(i,k)
  r = max(i,k)
  p = p+ishft(r*r-r,-1)
  q = min(j,l)
  s = max(j,l)
  q = q+ishft(s*s-s,-1)
  i1 = min(p,q)
  i2 = max(p,q)
  i1 = i1+ishft(i2*i2-i2,-1)
end

The Data Structure

Now that it is possible from (ij|kl) to obtain a unique 64-bit integer and the value of the integral, these can be used to store the data in a data structure similar to a hash table. In what follows, we call this data structure a map.

A map type is essentially an array of cache_maps. The cache_map type contains mainly a pointer to an array of values and a pointer to an array of keys, and a flag to know if the data is sorted by increasing keys.

The cache_map keys are stored on signed 16-bits integers, so the value of the key will not exceed 32768. This guarantees an upper bound on the total number of values of the cache_map, and the largest possible size of the keys array will be 64kiB. Most of the time, it is expected to fit in the L1 cache.

Search

When accessing the map data structure, the key is split in two. The 49 leading bits of the 64-bit key are used to extract the index of the cache map in the array of cache_maps. The 15 least significant bits are extracted to generate the cache_map key:

     000110101001001110011010010110100011101000011010  |  010011101001101
   <-------------- index of the cache_map ------------> <--- cache_key --->
                         integer*8                           integer*2

A direct access is made to the cache_map (usually in RAM, not in a cache). If the cache_map is not yet sorted, it is sorted. Then the cache_key is searched for in the array of cache_keys using a binary search. As the maximum possible size of the array is 32768, the data will be found in less than 16 tries, which are likely not to cross a memory page boundary. When the key is found, the corresponding value is at the same index in the array of values. If the key is not found, a value of zero is returned.

Insert

The key is appended to the array of keys and the value is appended to the array of values. When the cache_map arrays are completely filled with data, the size of the array is doubled (cache_map_reallocate), the keys are sorted (cache_map_sort), and the order is applied to the values array (dset_order). When duplicates are found (cache_map_unique), their values are added together such that there is always a unique key for a unique value. Upon addition of duplicates, the resulting value may be lower than a threshold so in that case the entry is removed from the cache_map (cache_map_shrink).

Thread safety

When the atomic integrals are written, they are not likely to access simultaneously the same cache_map. So the cache_maps are filled and sorted concurrently with almost no overhead. When all the data has been added, all the cache_maps are sorted concurrently.

During the 4-index transformation, the molecular integrals need to be frequently updated. Using this data structure, it is possible to modify the elements in a parallel section each cache_map has its own lock which is acquired when the cache_map is modified.