Como implementar a função de indexação eficiente para integrais de duas partículas <ij | kl>?

11

Este é um problema simples de enumeração de simetria. Eu dou o histórico completo aqui, mas nenhum conhecimento de química quântica é necessário.

O integrante dois partícula é: i j | k l = ψ * i ( x ) ψ * j ( x ' ) ψ k ( x ) ψ l ( x ' )ij|kl E tem os seguintes 4 simetrias: i j | k l = j i | l k = k l | i j = l k | j i I têm uma função que calcula o integral e armazena-os em uma matriz 1D, indexado como se segue:

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

onde a função ijkl2intindex2retorna um índice exclusivo, levando em consideração as simetrias acima. O único requisito é que, se você percorrer todas as combinações de i, j, k, l (de 1 a n cada), ele preencherá a int2matriz consecutivamente e atribuirá o mesmo índice a todas as combinações de ijkl relacionadas pelas opções acima. 4 simetrias.

Minha implementação atual no Fortran está aqui . Está muito lento. Alguém sabe como fazer isso de forma eficaz? (Em qualquer idioma.)

ψ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
fonte
d3x
1
d3xdx1dx2dx3x=(x1,x2,x3)d3x
xx=(x1,x2,x3)dx
d3x

Respostas:

5

[Editar: a quarta vez é o charme, finalmente algo sensato]

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

A combinação de ambos fornece o conjunto completo; portanto, juntar os dois loops nos fornece o conjunto completo de índices.

n

Em python, podemos escrever o seguinte iterador para fornecer os valores idx e i, j, k, l para cada cenário 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)

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

E então faça um loop sobre ele assim:

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
fonte
Olá Phil, muito obrigado pela resposta! Eu testei e há dois problemas. Por exemplo, idx_all (1, 2, 3, 4, 4) == idx_all (1, 2, 4, 3, 4) = 76. Mas <12 | 34> / = <12 | 43>. Só é igual se os orbitais forem reais. Portanto, sua solução parece ser o caso de 8 simetrias (veja meu exemplo do Fortran acima para uma versão mais simples, o ijkl2intindex ()). O segundo problema é que os índices não são consecutivos, colei os resultados aqui: gist.github.com/2703756 . Aqui estão os resultados corretos da minha rotina ijkl2intindex2 () acima: gist.github.com/2703767 .
Ondřej Čertík
1
@ OndřejČertík: Você quer um sinal associado? Faça com que o idxpair retorne um sinal se você trocou de ordem.
Deathbreath
OndřejČertík: Eu vejo a diferença agora. Como aponta @Deathbreath, você pode negar o índice, mas isso não será tão limpo para o loop geral. Vou pensar e atualizá-lo.
Phil H
Na verdade, negar o índice não funcionará totalmente, pois o idxpair errará o valor.
Phil H
<ij|kl>=<ji|kl>=<ij|lk>=<ji|lk>
ijkl[idxpair(indexij,indexkl,,)]signijsignkl
3

Aqui está uma idéia do uso de uma curva simples de preenchimento de espaço modificada para retornar a mesma chave para os casos de simetria (todos os trechos de código estão em 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:

  • O exemplo é python, mas se você incorporar as funções em seu código fortran e desenrolar seu loop interno para (i, j, k, l), deverá obter um desempenho decente.
  • Você pode calcular a chave usando flutuadores e, em seguida, converter a chave em número inteiro para usar como índice, isso permitiria ao compilador usar as unidades de ponto flutuante (por exemplo, o AVX está disponível).
  • Se N é uma potência de 2, então as multiplicações seriam apenas pequenas mudanças.
  • O tratamento para as simetrias não é eficiente na memória (ou seja, não produz uma indexação contínua) e utiliza cerca de 1/4 das entradas totais da matriz de índices.

Aqui está um exemplo de teste 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

Saída 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

Se for de interesse, a função inversa de forge_key é:

# 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
fonte
Você quis dizer "se n é uma potência de 2" em vez de múltiplos de 2?
Aron Ahmadia 17/05
sim, obrigado Aron. Eu escrevi essa resposta antes de ir jantar e Hulk estava escrevendo.
Fcruz 17/05/12
Esperto! No entanto, o índice máximo não é n ^ 4 (ou n ^ 4-1 se você começar de 0)? O problema é que, para o tamanho básico que eu quero poder, não cabe na memória. Com índice consecutivo, o tamanho da matriz é n ^ 2 * (n ^ 2 + 3) / 4. Hm, isso significa apenas 1/4 do tamanho total. Então talvez eu não deva me preocupar com o fator 4 no consumo de memória. Ainda assim, deve haver alguma maneira de codificar o índice consecutivo correto usando apenas essas 4 simetrias (melhor que minha solução feia no meu post, onde eu preciso fazer loops duplos).
Ondřej Čertík
Sim, está certo! Não sei como resolver com elegância (sem classificar e renumerar) o índice, mas o termo principal no uso de memória é O (N ^ 4). O fator de 4 deve fazer uma pequena diferença na memória para grande N.
fcruz
0

Isso não é apenas a generalização do problema de indexação de matriz simétrica compactada? A solução aí é offset (i, j) = i * (i + 1) / 2 + j, não é? Você não pode dobrar isso e indexar uma matriz 4D duplamente simétrica? A implementação que requer ramificação parece desnecessária.

Jeff
fonte