Как реализовать эффективную функцию индексации для двухчастичных интегралов <ij | kl>?

11

Это простая задача перечисления симметрии. Я даю полный обзор, но знание квантовой химии не требуется.

Двухчастичный интеграл : И имеет следующие 4 симметрии: У меня есть функция, которая вычисляет интеграл и сохраняет его в одномерном массиве , проиндексированном следующим образом:я J | к л = г | * я ( х ) г | * J ( х ' ) г | к ( х ) г | л ( х ' )ij|klя J | к л = J я | л к = к л | я J = л к | J я

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

где функция ijkl2intindex2возвращает уникальный индекс, принимая во внимание вышеуказанные симметрии. Единственное требование состоит в том, что если вы перебираете все комбинации i, j, k, l (от 1 до n каждая), он будет int2последовательно заполнять массив и назначать один и тот же индекс всем комбинациям ijkl, которые связаны с вышеуказанным 4 симметрии.

Моя текущая реализация в Фортране здесь . Это очень медленно. Кто-нибудь знает, как сделать это эффективно? (На любом языке.)

Подсказка: если орбитали действительны, то в дополнение к вышеуказанным симметриям можно поменять местами и чтобы мы получили всего 8 симметрий: и тогда можно реализовать очень быструю функцию для ее индексации, см. мою реализацию здесь . Я хотел бы найти эффективную схему индексации для случаев, когда орбитали не являются реальными.я к J л я J | к л = J яψi(x)ikjl

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

Примечание: функции, которые я реализовал, фактически принимают четыре числа , , , в так называемой «химической» записи , то есть аргументы и взаимозаменяемы, но это не важно.J K L ( я J | к л ) = я K | J л J Kijkl(ij|kl)=ik|jljk

Ондржей Чертик
источник
d3x
1
d3xdx1dx2dx3x=(x1,x2,x3)d3x
xx=(x1,x2,x3)dx
d3x

Ответы:

5

[Редактировать: 4-й раз это очарование, наконец-то что-то разумное]

nn2(n2+3)t(t(n))+t(t(n1))t(a)a, . Получив это, мы должны знать почему.t(a)=a(a+1)/2

Первый член проще - пары пар, где ij и , где - треугольный индекс . Это выполняется функцией, подобной этой:t i d ( a , b ) a , btid(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

Комбинация обоих из них дает полный набор, поэтому объединение обоих циклов дает нам полный набор индексов.

n

В python мы можем написать следующий итератор, чтобы дать нам значения idx и i, j, k, l для каждого отдельного сценария:

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

И затем зациклите это так:

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
Фил Х
источник
Привет Фил, большое спасибо за ответ! Я проверил это, и есть две проблемы. Например, idx_all (1, 2, 3, 4, 4) == idx_all (1, 2, 4, 3, 4) = 76. Но <12 | 34> / = <12 | 43>. Равен только если орбитали реальны. Таким образом, ваше решение похоже на случай 8 симметрий (см. Мой пример на Fortran выше для более простой версии, ijkl2intindex ()). Вторая проблема заключается в том, что индексы не являются последовательными, я вставил результаты здесь: gist.github.com/2703756 . Вот правильные результаты из моей процедуры ijkl2intindex2 (), приведенной выше: gist.github.com/2703767 .
Ондржей Чертик
1
@ OndřejČertík: Вы хотите знак, связанный? пусть idxpair вернет знак, если вы поменяли порядок.
Смертельное дыхание
Ондржейчертик: Теперь я вижу разницу. Как указывает @Deathbreath, вы можете отменить индекс, но это не будет таким же чистым для всего цикла. Я подумаю и обновлю это.
Фил Х
На самом деле отрицание индекса не будет работать полностью, так как idxpair получит неверное значение.
Фил Х
<ij|kl>=<ji|kl>=<ij|lk>=<ji|lk>
ijkl[idxpair(indexij,indexkl,,)]signijsignkl
3

Вот идея использования простой кривой заполнения пространства, модифицированной так, чтобы она возвращала один и тот же ключ для случаев симметрии (все фрагменты кода в 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)) 

Заметки:

  • Примером является python, но если вы встроите функции в код Fortran и развернете свой внутренний цикл для (i, j, k, l), вы должны получить приличную производительность.
  • Вы можете вычислить ключ, используя числа с плавающей запятой, а затем преобразовать ключ в целое число для использования в качестве индекса, это позволит компилятору использовать единицы с плавающей запятой (например, доступен AVX).
  • Если N - степень 2, то умножения будут просто сдвигами битов.
  • Обработка симметрий неэффективна в памяти (т. Е. Она не производит непрерывную индексацию) и использует около 1/4 всех записей массива индекса.

Вот тестовый пример для 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

Выход для 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

Если интересно, обратная функция 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
источник
Вы имели в виду «если n - степень 2» вместо кратного 2?
Арон Ахмадиа
да, спасибо Арон. Я написал этот ответ как раз перед обедом, и Халк писал.
fcruz
Умная! Однако не является ли максимальный индекс n ^ 4 (или n ^ 4-1, если вы начинаете с 0)? Проблема в том, что для базового размера, который я хочу иметь, он не помещается в память. При последовательном индексе размер массива равен n ^ 2 * (n ^ 2 + 3) / 4. Хм, в любом случае это всего лишь около 1/4 от полного размера. Так что, возможно, мне не следует беспокоиться о факторе 4 в потреблении памяти. Тем не менее, должен быть какой-то способ кодировать правильный последовательный индекс, используя только эти 4 симметрии (лучше, чем мое уродливое решение в моем посте, где мне нужно сделать двойные циклы).
Ондржей Чертик
да, это правильно! Я не знаю, как элегантно решить (без сортировки и перенумерации) индекс, но главный термин в использовании памяти - O (N ^ 4). Фактор 4 должен иметь небольшое различие в памяти для большого N.
fcruz
0

Разве это не просто обобщение задачи индексации упакованной симметричной матрицы? Там есть смещение (i, j) = i * (i + 1) / 2 + j, не так ли? Разве вы не можете удвоить это и проиндексировать двусимметричный 4D массив? Реализация, требующая ветвления, кажется ненужной.

Джефф
источник