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

11

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 ' )ij|kl 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:

ij|kl=ψi(x)ψj(x)ψk(x)ψl(x)|xx|d3xd3x
ij|kl=ji|lk=kl|ij=lk|ji
int2
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.)

ψi(x)ikjl

ij|kl=ji|lk=kj|il=il|kj=
=kl|ij=lk|ji=il|kj=kj|il

ijkl(ij|kl)=ik|jljk

Ondřej Čertík
la source
d3x
1
d3xdx1dx2dx3x=(x1,x2,x3)d3x
xx=(x1,x2,x3)dx
d3x

Réponses:

5

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

nn2(n2+3)t(t(n))+t(t(n1))t(a)at(a)=a(a+1)/2

ijtid(i,j)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

llk

t(t(n1))

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)

in3+jn2+kn+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
Phil H
la source
Salut Phil, merci beaucoup pour la réponse! Je l'ai testé et il y a deux problèmes. Par exemple idx_all (1, 2, 3, 4, 4) == idx_all (1, 2, 4, 3, 4) = 76. Mais <12 | 34> / = <12 | 43>. Elle n'est égale que si les orbitales sont réelles. Votre solution semble donc être pour le cas de 8 symétries (voir mon exemple Fortran ci-dessus pour une version plus simple, ijkl2intindex ()). Le deuxième problème est que les indices ne sont pas consécutifs, j'ai collé les résultats ici: gist.github.com/2703756 . Voici les résultats corrects de ma routine ijkl2intindex2 () ci-dessus: gist.github.com/2703767 .
Ondřej Čertík
1
@ OndřejČertík: Vous voulez un signe associé? demandez à idxpair de retourner un signe si vous avez changé de commande.
Deathbreath
OndřejČertík: Je vois maintenant la différence. Comme le souligne @Deathbreath, vous pouvez annuler l'index, mais ce ne sera pas aussi propre pour la boucle globale. Je vais réfléchir et le mettre à jour.
Phil H
En fait, la négation de l'index ne fonctionnera pas complètement car idxpair obtiendra une valeur incorrecte.
Phil H
@PhilH: Ne niez pas l'index j
<ij|kl>=<ji|kl>=<ij|lk>=<ji|lk>
ijkl[idxpair(indexij,indexkl,,)]signijsignkl
3

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)
fcruz
la source
Voulez-vous dire "si n est une puissance de 2" au lieu de multiples de 2?
Aron Ahmadia
oui, merci Aron. J'ai écrit cette réponse juste avant d'aller dîner et Hulk écrivait.
fcruz
Intelligent! Cependant, l'indice maximal n'est-il pas n ^ 4 (ou n ^ 4-1 si vous partez de 0)? Le problème est que pour la taille de base que je veux pouvoir faire, elle ne rentrera pas dans la mémoire. Avec un index consécutif, la taille du tableau est n ^ 2 * (n ^ 2 + 3) / 4. Hm, ce n'est de toute façon qu'environ 1/4 de la taille complète. Alors peut-être que je ne devrais pas m'inquiéter du facteur 4 de la consommation de mémoire. Cependant, il doit y avoir un moyen de coder l'index consécutif correct en utilisant uniquement ces 4 symétries (mieux que ma vilaine solution dans mon article, où je dois faire des boucles doubles).
Ondřej Čertík
oui, c'est vrai! Je ne sais pas comment résoudre élégamment (sans trier ni renuméroter) l'index mais le terme principal dans l'utilisation de la mémoire est O (N ^ 4). Le facteur 4 devrait faire une petite différence dans la mémoire pour les grands N.
fcruz
0

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.

Jeff
la source