计算科学 - 如何实现两个粒子积分<ij|kl>的高效索引函数？ - 吾爱随笔录

如何实现两个粒子积分<ij|kl>的高效索引函数？

计算科学算法表现

2021-12-07 04:19:07

这是一个简单的对称枚举问题。我在这里给出了完整的背景，但不需要量子化学知识。

两粒子积分 $\langle ij|kl\rangle$ 是：

⟨ 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'$ 它具有以下 4 个对称性：

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

其中函数ijkl2intindex2返回一个唯一索引，考虑到上述对称性。唯一的要求是，如果您遍历 i、j、k、l 的所有组合（每个从 1 到 n），它将int2连续填充数组，并将相同的索引分配给与上述相关的所有 ijkl 组合4 对称性。

我目前在 Fortran 中的实现在这里。它非常慢。有人知道如何有效地做到这一点吗？（任何语言。）

提示：如果轨道 $\psi_i({\bf x})$ 是真实的，那么除了上述对称性之外，还可以交换 $i\leftrightarrow k$ 和 $j\leftrightarrow l$ 所以我们总共得到8个对称性：

⟨ 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$ 然后可以实现一个非常快速的函数来索引它，请参阅我的实现here。我想为轨道不真实的情况找到一些有效的索引方案。

注意：我实现的函数实际上接受四个数字 $i$ , $j$ , $k$ , $l$ 在所谓的“化学”符号中 $(ij|kl) = \langle ik|jl\rangle$ ，即 $j$ 和 $k$ 参数可以互换，但这并不重要。

3个回答

[编辑：第四次的魅力，终于有点明智]

我倒退了：我从另一个显示基于过滤器的方法的答案开始，并用它来生成一系列值的所有有效组合 $n$ ，并在在线整数序列数据库中查找序列。组合数为 $n^2(n^2+3)$ ，这看起来不太可能（为什么是 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$ 没有 if/skip 的循环来检查三角形索引。

第二学期 $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

这两者的组合给出了完整的集合，因此将两个循环放在一起给出了完整的索引集合。

一个重要的问题是，对于任意 i,j,k,l，这些模式很难计算。所以我会建议一张地图，它产生给定 i,j,k,l 的索引。坦率地说，如果您要这样做，您不妨使用 generate+filter 方法，因为您只需要对给定的情况执行一次 $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)

在 fortran 中，我们只需要运行循环并存储值。我们可以使用一个简单的索引将 i,j,k,l 组合存储为单个值 ( $i n^3 + j n^2 + k n + l$ )，并将这些值存储在索引与上述索引相同的数组中。然后我们可以遍历这个数组并从值中检索 i,j,k,l。要获得任意 i,j,k,l 的 idx，需要一个反向映射和一个过滤器来处理对称性，尽管我们可能可以从上述结构构造一个函数。fortran 中的 idx 数组生成函数为：

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

这是使用修改后的简单空间填充曲线以针对对称情况返回相同键的想法（所有代码片段都在 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)

这不仅仅是压缩对称矩阵索引问题的推广吗？那里的解决方案是 offset(i,j) = i*(i+1)/2+j，不是吗？你不能加倍努力并索引一个双对称 4D 数组吗？需要分支的实现似乎没有必要。

其它你可能感兴趣的问题

上一篇BFGS 对初始 Hessian 近似的敏感性下一篇相对于六面体单元的非结构化网格对点云进行排序