Jak wdrożyć wydajną funkcję indeksowania dla dwóch całek cząstek <ij | kl>?

11

Jest to prosty problem z wyliczaniem symetrii. Podaję tutaj pełne tło, ale nie jest wymagana znajomość chemii kwantowej.

Integralną dwóch cząstek jest: ja j | K L = ψ * I ( x ) ψ * j ( x ' ) ψ k ( x ) ψ L ( x ' )ij|kl I ma następujące 4 symetrie: í j | k l = j I | l k = k l | I j = l k | j I Mam funkcję, która oblicza całki i zapisuje je w tablicy 1D, indeksowane, co następuje:

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

gdzie funkcja ijkl2intindex2zwraca unikalny indeks, biorąc pod uwagę powyższe symetrie. Jedynym wymaganiem jest to, że jeśli zapętlisz wszystkie kombinacje i, j, k, l (od 1 do n każda), wypełni ona int2tablicę kolejno i przypisze ten sam indeks do wszystkich kombinacji ijkl, które są powiązane powyższym 4 symetrie.

Moja obecna implementacja w Fortran jest tutaj . To jest bardzo wolne. Czy ktoś wie, jak to zrobić skutecznie? (W dowolnym języku.)

ψ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
źródło
d3x
1
d3xdx1dx2dx3x=(x1,x2,x3)d3x
xx=(x1,x2,x3)dx
d3x

Odpowiedzi:

5

[Edytuj: Czwarty raz to urok, w końcu coś sensownego]

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

Połączenie obu tych elementów daje pełny zestaw, więc połączenie obu pętli daje nam pełny zestaw wskaźników.

n

W pythonie możemy napisać następujący iterator, aby podać nam wartości idx i i, j, k, l dla każdego innego scenariusza:

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

A następnie zapętlić go tak:

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.
źródło
Cześć Phil, wielkie dzięki za odpowiedź! Przetestowałem to i są dwa problemy. Na przykład idx_all (1, 2, 3, 4, 4) == idx_all (1, 2, 4, 3, 4) = 76. Ale <12 | 34> / = <12 | 43>. Jest równy tylko wtedy, gdy orbitale są prawdziwe. Więc twoje rozwiązanie wydaje się być w przypadku 8 symetrii (zobacz mój przykład Fortran powyżej, aby uzyskać prostszą wersję, ijkl2intindex ()). Drugi problem polega na tym, że wskaźniki nie są kolejne, wkleiłem wyniki tutaj: gist.github.com/2703756 . Oto poprawne wyniki z powyższej procedury ijkl2intindex2 (): gist.github.com/2703767 .
Ondřej Čertík
1
@ OndřejČertík: Chcesz skojarzyć znak? niech idxpair zwróci znak, jeśli zmieniłeś zamówienie.
Deathbreath
OndřejČertík: Teraz widzę różnicę. Jak wskazuje @Deathbreath, możesz zanegować indeks, ale nie będzie to tak czyste dla całej pętli. Zastanowię się i zaktualizuję.
Phil H
W rzeczywistości zanegowanie indeksu nie będzie w pełni działać, ponieważ idxpair błędnie poda wartość.
Phil H
<ij|kl>=<ji|kl>=<ij|lk>=<ji|lk>
ijkl[idxpair(indexij,indexkl,,)]signijsignkl
3

Oto pomysł użycia prostej krzywej wypełniania spacji zmodyfikowanej w celu zwrócenia tego samego klucza dla przypadków symetrii (wszystkie fragmenty kodu znajdują się w pythonie).

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

Uwagi:

  • Przykładem jest python, ale jeśli wstawisz funkcje do kodu fortran i rozwiniesz wewnętrzną pętlę dla (i, j, k, l), powinieneś uzyskać przyzwoitą wydajność.
  • Możesz obliczyć klucz za pomocą liczb zmiennoprzecinkowych, a następnie przekonwertować klucz na liczbę całkowitą, aby użyć go jako indeksu, co pozwoliłoby kompilatorowi na użycie jednostek zmiennoprzecinkowych (np. AVX jest dostępny).
  • Jeśli N jest potęgą 2, mnożenia byłyby tylko przesunięciami bitowymi.
  • Traktowanie symetrii nie jest wydajne w pamięci (tj. Nie powoduje ciągłego indeksowania) i zużywa około 1/4 wszystkich wpisów tablicy indeksów.

Oto przykładowy test dla 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

Dane wyjściowe dla 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

Jeśli jest to interesujące, funkcją odwrotną do forge_key jest:

# 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
źródło
Czy chodziło Ci o „jeśli n jest potęgą 2” zamiast wielokrotności 2?
Aron Ahmadia
tak, dzięki Aron. Napisałem tę odpowiedź tuż przed pójściem na obiad, a Hulk pisał.
fcruz
Sprytny! Czy jednak nie jest maksymalny indeks n ^ 4 (lub n ^ 4-1, jeśli zaczynasz od 0)? Problem polega na tym, że dla rozmiaru podstawowego, który chcę mieć, nie zmieści się on w pamięci. Przy kolejnym indeksie rozmiar tablicy wynosi n ^ 2 * (n ^ 2 + 3) / 4. Hm, to i tak tylko około 1/4 pełnego rozmiaru. Więc może nie powinienem martwić się współczynnikiem 4 zużycia pamięci. Nadal jednak musi być jakiś sposób na zakodowanie poprawnego kolejnego indeksu przy użyciu tylko tych 4 symetrii (lepiej niż moje brzydkie rozwiązanie w moim poście, gdzie muszę robić podwójne pętle).
Ondřej Čertík
tak, właśnie tak! Nie wiem, jak elegancko rozwiązać (bez sortowania i numerowania) indeks, ale wiodącym terminem użycia pamięci jest O (N ^ 4). Współczynnik 4 powinien mieć niewielką różnicę w pamięci dla dużej N.
fcruz
0

Czy to nie tylko uogólnienie problemu upakowanego indeksowania macierzy symetrycznej? Rozwiązaniem jest przesunięcie (i, j) = i * (i + 1) / 2 + j, prawda? Nie możesz podwoić tego i zindeksować podwójnie symetryczną macierz 4D? Implementacja wymagająca rozgałęzień wydaje się niepotrzebna.

Jeff
źródło