Comment implémenter une fonction d'indexation efficace pour deux intégrales de particules <ij | kl>?

Il s'agit d'un simple problème d'énumération de symétrie. Je donne ici le contexte complet, mais aucune connaissance de la chimie quantique n'est nécessaire.

L'intégrale de deux particules est: ⟨ i j | k l ⟩ = ∫ ψ * i ( x ) ψ * j ( x ' ) ψ k ( x ) ψ l ( x '$\langle ij|kl\rangle$ Et il a les 4 symétries suivantes: ⟨ i j | k l ⟩ = ⟨ j i | l k ⟩ = ⟨ k l | i j ⟩ = ⟨ l k | j i ⟩ I ont une fonction qui calcule l'intégrale et les stocke dans un tableau 1D, indexées comme suit:

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

où la fonction ijkl2intindex2renvoie un index unique, en tenant compte des symétries ci-dessus. La seule exigence est que si vous bouclez sur toutes les combinaisons de i, j, k, l (de 1 à n chacune), il remplira le int2tableau consécutivement et attribuera le même index à toutes les combinaisons ijkl qui sont liées par ce qui précède 4 symétries.

Mon implémentation actuelle à Fortran est ici . C'est très lent. Quelqu'un sait-il comment le faire efficacement? (Dans n'importe quelle langue.)

$\psi_i({\bf x})$$i\leftrightarrow k$$j\leftrightarrow l$

$i$$j$$k$$l$$(ij|kl) = \langle ik|jl\rangle$$j$$k$

[Edit: la 4ème fois c'est le charme, enfin quelque chose de sensé]

$n$$n^2(n^2+3)$$t(t(n))+t(t(n-1))$$t(a)$$a$$t(a)=a(a+1)/2$

$i \ge j$$tid(i,j) \ge tid(k,l)$$tid(a,b)$$a,b$

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

$l$$l$$k$

$t(t(n-1))$

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

La combinaison des deux donne l'ensemble complet, donc mettre les deux boucles ensemble nous donne l'ensemble complet des indices.

$n$

En python, nous pouvons écrire l'itérateur suivant pour nous donner les valeurs idx et i, j, k, l pour chaque scénario différent:

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)

$i n^3 + j n^2 + k n + l$

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

Et puis bouclez-le ainsi:

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

Voici une idée d'utiliser une simple courbe de remplissage d'espace modifiée pour retourner la même clé pour les cas de symétrie (tous les extraits de code sont en 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)) 

Remarques:

• L'exemple est python mais si vous insérez les fonctions dans votre code fortran et déroulez votre boucle interne pour (i, j, k, l), vous devriez obtenir des performances décentes.
• Vous pouvez calculer la clé en utilisant des flottants puis convertir la clé en entier pour l'utiliser comme index, cela permettrait au compilateur d'utiliser les unités à virgule flottante (par exemple AVX est disponible).
• Si N est une puissance de 2, les multiplications ne seraient que des décalages de bits.
• Le traitement des symétries n'est pas efficace en mémoire (c'est-à-dire qu'il ne produit pas d'indexation continue) et utilise environ 1/4 du total des entrées du tableau d'index.

Voici un exemple de test pour 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

Sortie pour 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

Si cela vous intéresse, la fonction inverse de forge_key est:

# 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)

N'est-ce pas simplement la généralisation du problème d'indexation de matrice symétrique compactée? La solution y est offset (i, j) = i * (i + 1) / 2 + j, n'est-ce pas? Vous ne pouvez pas doubler cela et indexer un tableau 4D doublement symétrique? L'implémentation nécessitant une ramification semble inutile.