¿Cómo implementar una función de indexación eficiente para dos integrales de partículas <ij | kl>?

Este es un simple problema de enumeración de simetría. Doy todos los antecedentes aquí, pero no se necesita conocimiento de química cuántica.

La integral de dos partículas es: Y tiene las siguientes 4 simetrías: Tengo una función que calcula la integral y la almacena en una matriz 1D , indexada de la siguiente manera:⟨ i j | k l ⟩ = ∫ ψ * i ( x ) ψ * j ( x ' ) ψ k ( x ) ψ l ( x '$\langle ij|kl\rangle$ ⟨ i j | k l ⟩ = ⟨ j i | l k ⟩ = ⟨ k l | i j ⟩ = ⟨ l k | j i

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

donde la función ijkl2intindex2devuelve un índice único, teniendo en cuenta las simetrías anteriores. El único requisito es que si realiza un bucle sobre todas las combinaciones de i, j, k, l (de 1 a n cada una), llenará la int2matriz consecutivamente y asignará el mismo índice a todas las combinaciones de ijkl que están relacionadas por lo anterior 4 simetrías

Mi implementación actual en Fortran está aquí . Es muy lento. ¿Alguien sabe cómo hacer esto de manera efectiva? (En cualquier idioma.)

Sugerencia: si los orbitales son reales, además de las simetrías anteriores, se puede intercambiar y para obtener 8 simetrías en total: y luego uno puede implementar una función muy rápida para indexarlo, vea mi implementación aquí . Me gustaría encontrar un esquema de indexación eficiente para los casos en que los orbitales no son reales.$\psi_i({\bf x})$$i\leftrightarrow k$$j\leftrightarrow l$

Nota: las funciones que implementé en realidad aceptan los cuatro números , , , en una notación de "química" , es decir, los argumentos y se intercambian, Pero esto no es importante.j k l ( i j | k l ) = ⟨ i k | j l ⟩ j $i$$j$$k$$l$$(ij|kl) = \langle ik|jl\rangle$$j$$k$

[Editar: cuarta vez es el encanto, por fin algo sensato]

Llegué a esto al revés: comencé con otra respuesta que mostraba el enfoque basado en filtros, y usé eso para generar todas las combinaciones válidas para una serie de valores de , y busqué la secuencia en la base de datos de secuencias de enteros en línea. El número de combinaciones es , lo que parece poco probable (¿por qué el 3?). También es donde es el número triangular de , . Habiendo entendido esto, necesitamos saber por qué.$n$$n^2(n^2+3)$$t(t(n))+t(t(n-1))$$t(a)$$a$$t(a)=a(a+1)/2$

El primer término es más simple: los pares de pares donde y , donde es el índice triangular de . Eso se cumple con una función como esta:t i d ( i , j ) ≥ t i d ( k , l ) t i d ( a , b ) a , $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

donde el segundo bucle es porque no podemos anidar el bucle completo dentro del bucle sin un if / skip para verificar los índices triangulares.l $l$$l$$k$

El segundo término es donde el primer par es ascendente y el segundo par es descendente (nota, no ==, como se manejan arriba).$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 combinación de ambos proporciona el conjunto completo, por lo que al unir ambos bucles se obtiene el conjunto completo de índices.

Un problema importante es que estos patrones son difíciles de calcular para un arbitrario i, j, k, l. Por lo tanto, sugeriría un mapa que arroje el índice dado i, j, k, l. Francamente, si está haciendo esto, también podría usar el enfoque de generar + filtro, porque solo necesita hacerlo una vez para un dado . El lado positivo del método anterior es que al menos tiene una estructura de bucle predecible.$n$

En python podemos escribir el siguiente iterador para darnos los valores idx e i, j, k, l para cada escenario diferente:

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)

En fortran solo tendríamos que ejecutar el ciclo y almacenar los valores. Podemos usar un índice simple para almacenar la combinación i, j, k, l como un valor único ( ), y almacenar estos valores en una matriz cuyo índice es el mismo que el índice encima. Luego podemos iterar sobre esta matriz y recuperar i, j, k, l de los valores. Para obtener el idx para i, j, k, l arbitrario, se requeriría un mapa inverso y un filtro para manejar la simetría, aunque probablemente podríamos construir una función a partir de la estructura anterior. La función de generación de matriz idx en fortran sería:$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

Y luego repítelo así:

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

Aquí hay una idea de usar una curva de relleno de espacio simple modificada para devolver la misma clave para los casos de simetría (todos los fragmentos de código están 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)) 

Notas:

• El ejemplo es python, pero si integra las funciones en su código fortran y desenrolla su bucle interno para (i, j, k, l), debería obtener un rendimiento decente.
• Podría calcular la clave usando flotantes y luego convertir la clave en un entero para usarla como índice, esto permitiría al compilador usar las unidades de coma flotante (por ejemplo, AVX está disponible).
• Si N es una potencia de 2, entonces las multiplicaciones serían solo cambios de bits.
• El tratamiento de las simetrías no es eficiente en la memoria (es decir, no produce una indexación continua) y utiliza alrededor de 1/4 de las entradas de la matriz de índice total.

Aquí hay un ejemplo de prueba para 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

Salida para 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 es de interés, la función inversa de forge_key es:

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

¿No se trata solo de la generalización del problema de indexación de matriz simétrica empaquetada? La solución es offset (i, j) = i * (i + 1) / 2 + j, ¿no es así? ¿No puede duplicar esto e indexar una matriz 4D doblemente simétrica? La implementación que requiere ramificación parece innecesaria.