晶胞中的距离矩阵（考虑对称性）

Posted 2023-02-23

【中文标题】晶胞中的距离矩阵（考虑对称性）

【英文标题】：Distance matrix in unit cell (accounting for symmetry)

【发布时间】：2022-01-07 02:41:24

【问题描述】：

我在计算大距离矩阵时遇到了问题。然而，这是一个特定的距离矩阵：它是一个单元格中的点矩阵。该函数获取分数坐标（在所有维度上都在 0 和 1 之间），我想计算距离矩阵，因为在单位单元的每个邻居中都有一个相同的点副本，因此正确的距离可能与副本而不是与单元格内的其他点一起。

您知道是否可以为此使用 scipy 或 numpy 预编码的 C 库来做任何事情？我已经完成了一个 numba 代码，它可以工作但运行速度很慢。这里我有一个 13160 个点的列表，我想计算一个 13160*13160 的距离矩阵，即包含 173185600 个元素。

原理是：对于每个坐标，计算第一个点与第二个点的平方分数距离，或者在单元格内，或者在它的两个邻居之一（后面和前面）。然后得到每个坐标的平方距离的最小值，得到与笛卡尔坐标对应的欧几里得距离。

当前耗时为：40.82661843299866秒

你知道我是否可以通过任何方式让它运行得更快，还是我的数据集太大而没有更多的事情要做？

下面是代码：

def getDistInCell(fract, xyz, n_sg, a, b, c):                           #calculates the distance matrix accounting for symmetry   
    dist = np.zeros((n_sg, n_sg))
    for i in range(n_sg):
        for j in range(n_sg):
            #we evaluate the closest segment according to translation to neighbouring cells
            diff_x = np.zeros((3))
            diff_y = np.zeros((3))
            diff_z = np.zeros((3))
            
            diff_x[0] = (fract[i][0] - (fract[j][0] - 1))**2
            diff_x[1] = (fract[i][0] - (fract[j][0]    ))**2
            diff_x[2] = (fract[i][0] - (fract[j][0] + 1))**2
            
            diff_y[0] = (fract[i][1] - (fract[j][1] - 1))**2
            diff_y[1] = (fract[i][1] - (fract[j][1]    ))**2
            diff_y[2] = (fract[i][1] - (fract[j][1] + 1))**2
            
            diff_z[0] = (fract[i][2] - (fract[j][2] - 1))**2
            diff_z[1] = (fract[i][2] - (fract[j][2]    ))**2
            diff_z[2] = (fract[i][2] - (fract[j][2] + 1))**2
            
            #get correct shifts
            shx = np.argmin(diff_x) - 1
            shy = np.argmin(diff_y) - 1
            shz = np.argmin(diff_z) - 1
            
            #compute cartesian distance
            dist[i][j] = np.sqrt((xyz[i][0] - (xyz[j][0] + shx * a)) ** 2 + (xyz[i][1] - (xyz[j][1] + shy * b)) ** 2 + (xyz[i][2] - (xyz[j][2] + shz * c)) ** 2)
    
    return dist

【问题讨论】：

请为代码提供一些（确定性的）输入。此外，最好有一个明确的dist[i][j] 公式（而不是试图从代码中提取它）。有数学解释吗？

***.com/questions/52649815/… 和 ***.com/questions/67035793/…

【参考方案1】：

这是一个基于 BallTree 的解决方案草图

我创建随机点，13160

import numpy as np

n=13160 
np.random.seed(1)

points=np.random.uniform(size=(n,3))

创建镜像/对称，例如

from itertools import product

def create_symmetries( points ):
    
    symmetries = []
    
    for sym in product([0,-1,1],[0,-1,1],[0,-1,1]):
        new_symmetry = points.copy()
        
        diff_x, diff_y, diff_z = sym
        new_symmetry[:,0] = new_symmetry[:,0] + diff_x
        new_symmetry[:,1] = new_symmetry[:,1] + diff_y
        new_symmetry[:,2] = new_symmetry[:,2] + diff_z

        symmetries.append(new_symmetry)
    
    return symmetries

并创建更大的数据集，包括对称性；

all_symmetries = np.concatenate( create_symmetries(points) )

要获得最接近的，请使用k=2，因为最接近的是点本身，第二接近的是最接近的对称（包括它自己的，所以要小心）

%%time 
import numpy as np
from sklearn.neighbors import BallTree


tree = BallTree(all_symmetries, leaf_size=15, metric='euclidean')

dist, idx = tree.query(points, k=2, return_distance=True)

这需要

CPU times: user 275 ms, sys: 2.77 ms, total: 277 ms
Wall time: 275 ms

【讨论】：

抱歉，我不明白 Neighbours_1[:,1] -1 究竟做了什么。你能帮我分解一下吗？我需要以下副本： - 左侧的单元格，右侧的单元格 - 顶部的单元格，底部的单元格 - 前面的单元格，背面的单元格

我已经更新了 symmetires - 但不确定这些是否全部。想法是您只需复制数据并检查最接近的...