Bagaimana menerapkan fungsi pengindeksan yang efisien untuk dua integral partikel <ij | kl>?

Ini adalah masalah enumerasi simetri sederhana. Saya memberikan latar belakang penuh di sini, tetapi tidak ada pengetahuan tentang kimia kuantum diperlukan.

Kedua partikel terpisahkan adalah: $\langle ij|kl\rangle$ Dan memiliki berikut 4 simetri: saya fungsi yang menghitung integral dan menyimpannya dalam array 1D, diindeks sebagai berikut:

⟨ i j | k l ⟩ = \int \frac{ψ_{i}^{*} (x) ψ_{j}^{*} (x^{'}) ψ_{k} (x) ψ_{l} (x^{'})}{| x - x^{'} |} d^{3} x d^{3} x^{'}

$\langle ij|kl\rangle = \int {\psi_i^*({\bf x})\psi_j^*({\bf x}')\psi_k({\bf x})\psi_l({\bf x}') \over | {\bf x} - {\bf x}' |}\, d^3 x \,d^3 x'$

⟨ i j | k l ⟩ = ⟨ j i | l k ⟩ = ⟨ k l | i j ⟩ = ⟨ l k | j i ⟩

$\langle ij|kl\rangle = \langle ji|lk\rangle = \langle kl|ij\rangle = \langle lk|ji\rangle$ int2

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

di mana fungsi ijkl2intindex2mengembalikan indeks unik, dengan mempertimbangkan simetri di atas. Satu-satunya persyaratan adalah bahwa jika Anda mengulang semua kombinasi i, j, k, l (masing-masing dari 1 hingga n), ia akan mengisi int2array secara berurutan dan akan menetapkan indeks yang sama untuk semua kombinasi ijkl yang terkait dengan hal di atas 4 simetri.

Implementasi saya saat ini di Fortran ada di sini . Sangat lambat. Adakah yang tahu cara melakukan ini secara efektif? (Dalam bahasa apa pun.)

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

⟨ i j | k l ⟩ = ⟨ j i | l k ⟩ = ⟨ k j | i l ⟩ = ⟨ i l | k j ⟩ =

$\langle ij|kl\rangle = \langle ji|lk\rangle = \langle kj|il\rangle = \langle il|kj\rangle =$

= ⟨ k l | i j ⟩ = ⟨ l k | j i ⟩ = ⟨ i l | k j ⟩ = ⟨ k j | i l ⟩

$= \langle kl|ij\rangle = \langle lk|ji\rangle = \langle il|kj\rangle = \langle kj|il\rangle$

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

algorithms performance Ondřej Čertík
sumber

d^{3} x

$d^3x$

d^{3} x

$d^3 x$

d x_{1} d x_{2} d x_{3}

$d x_1 d x_2 d x_3$

x = (x_{1}, x_{2}, x_{3})

${\bf x} = (x_1, x_2, x_3)$

d^{3} x

$d^3 x$

x

$x$

x = (x_{1}, x_{2}, x_{3})

$\mathbf{x} =(x_1,x_2,x_3)$

d x

$\mathrm{d}\mathbf{x}$

d^{3} x

$d^3 x$

Jawaban:

[Sunting: ke-4 kalinya pesona, akhirnya sesuatu yang masuk akal]

$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

Kombinasi keduanya memberikan set lengkap, jadi menempatkan kedua loop bersama memberi kita set lengkap indeks.

$n$

Dalam python kita dapat menulis iterator berikut untuk memberi kita nilai idx dan i, j, k, l untuk setiap skenario yang berbeda:

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

Dan kemudian lilitkan seperti itu:

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
sumber

Hai Phil, terima kasih banyak atas jawabannya! Saya mengujinya, dan ada dua masalah. Misalnya idx_all (1, 2, 3, 4, 4) == idx_all (1, 2, 4, 3, 4) = 76. Tetapi <12 | 34> / = <12 | 43>. Itu sama hanya jika orbital itu nyata. Jadi solusi Anda tampaknya untuk kasus 8 simetri (lihat contoh Fortran saya di atas untuk versi yang lebih sederhana, ijkl2intindex ()). Masalah kedua adalah bahwa indeks tidak berurutan, saya menempelkan hasilnya di sini: gist.github.com/2703756 . Berikut adalah hasil yang benar dari rutin ijkl2intindex2 () saya di atas: gist.github.com/2703767 .

OndČej Čertík

@ OndřejČertík: Anda ingin tanda dikaitkan? minta idxpair mengembalikan tanda jika Anda mengganti pesanan.

Deathbreath

OndřejČertík: Saya melihat perbedaannya sekarang. Seperti yang ditunjukkan oleh @Deathbreath, Anda dapat meniadakan indeks, tetapi itu tidak akan sebersih untuk keseluruhan loop. Saya akan berpikir dan memperbaruinya.

Phil H

Sebenarnya, meniadakan indeks tidak akan berfungsi sepenuhnya karena idxpair akan mendapatkan nilai yang salah.

Phil H

< i j | k l >= - < j i | k l >= - < i j | l k >=< j i | l k >

$<ij|kl>=-<ji|kl>=-<ij|lk> =<ji|lk>$

i j k l [i d x p a i r (i n d e x_{i j}, i n d e x_{k l},,)] s i g n_{i j} s i g n_{k l}

$ijkl[idxpair(index_{ij},index_{kl},,)]sign_{ij}sign_{kl}$

Berikut ini adalah ide untuk menggunakan kurva pengisian ruang sederhana yang dimodifikasi untuk mengembalikan kunci yang sama untuk kasus simetri (semua cuplikan kode menggunakan 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))

Catatan:

Contohnya adalah python tetapi jika Anda memasukkan fungsi ke dalam kode fortran Anda dan membuka gulungan lingkaran dalam untuk (i, j, k, l), Anda harus mendapatkan kinerja yang layak.
Anda dapat menghitung kunci menggunakan float dan kemudian mengubah kunci menjadi integer untuk digunakan sebagai indeks, ini akan memungkinkan kompiler untuk menggunakan unit floating point (mis. AVX tersedia).
Jika N adalah kekuatan 2 maka penggandaannya hanya akan sedikit bergeser.
Perawatan untuk simetri tidak efisien dalam memori (yaitu tidak menghasilkan pengindeksan berkelanjutan) dan menggunakan sekitar 1/4 dari total entri array indeks.

Ini adalah contoh uji untuk 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

Output untuk 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

Jika menarik, fungsi terbalik dari forge_key adalah:

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

Apakah maksud Anda "jika n adalah kekuatan 2", bukan kelipatan 2?

Aron Ahmadia

ya terima kasih Aron. Saya menulis jawaban ini sebelum pergi makan malam dan Hulk menulis.

fcruz

Pintar! Namun, bukankah indeks maksimum n ^ 4 (atau n ^ 4-1 jika Anda mulai dari 0)? Masalahnya adalah bahwa untuk ukuran dasar yang saya ingin dapat lakukan, itu tidak muat ke dalam memori. Dengan indeks berurutan, ukuran array adalah n ^ 2 * (n ^ 2 + 3) / 4. Hm, itu hanya sekitar 1/4 dari ukuran penuh. Jadi mungkin saya tidak perlu khawatir tentang faktor 4 dalam konsumsi memori. Meski begitu, harus ada beberapa cara untuk menyandikan indeks berurutan yang benar hanya menggunakan 4 simetri ini (lebih baik daripada solusi jelek saya di posting saya, di mana saya perlu melakukan loop ganda).

Ondřej Čertík

ya itu benar! Saya tidak tahu bagaimana menyelesaikan secara elegan (tanpa memilah dan memberi nomor baru) indeks tetapi istilah utama dalam penggunaan memori adalah O (N ^ 4). Faktor 4 harus membuat perbedaan kecil dalam memori untuk N. besar

fcruz

Apakah ini bukan hanya generalisasi dari masalah pengindeksan matriks simetris dikemas? Solusinya ada offset (i, j) = i * (i + 1) / 2 + j, bukan? Tidak bisakah Anda menggandakan ini dan mengindeks array 4D ganda simetris? Implementasi yang membutuhkan percabangan tampaknya tidak perlu.

Jeff
sumber