为数组的一个元素寻找邻居并确保它有 <4 个答案

【问题标题】：Finding neighbours for a an element of an array and making sure it has <4 of them为数组的一个元素寻找邻居并确保它有 <4 个
【发布时间】：2016-08-25 00:20:59
【问题描述】：

我有一个原子坐标数组（钻石晶格）看起来像这样：

                      [[0.0, 0.0, 0.0],
                      [0.0, 0.5, 0.5],
                      [0.5, 0.0, 0.5],
                      [0.5, 0.5, 0.0],
                      [0.25, 0.25, 0.25],
                      [0.25, 0.75, 0.75],
                      [0.75, 0.25, 0.75],
                      [0.75, 0.75, 0.25]]

每个原子都有一个 id。我需要创建一个键列表，这是一个给定原子所连接的原子的邻居列表。所有内部原子必须有 4 个邻居（键）。所有表面原子不能有超过 4 个邻居（键）。

我知道如何处理内部原子：使用快速且运行良好的 kdTree。

kdtree = scipy.spatial.KDTree(Atom_xyz[:,])
Neighs = kdtree.query_pairs(RCUT=1.0)
Neighs = list(Neighs)

但是它不适用于表面原子，因为它不适用于周期性图像，我需要确保我的邻居不超过 4 个。我想知道你们是否可以帮助我提供一个算法。

在这里我附上我的晶格的图像，其中原子被可视化为青色，而键则是弯曲的白色。

只有内部原子结合的金刚石晶格：

【问题讨论】：

标签： python numpy scipy simulation kdtree

【解决方案1】：

我刚带了四个邻居：

注意：内部原子可能存在Voronoi图没有四个邻居的病态情况）

（我用http://docs.enthought.com/mayavi/mayavi/auto/example_chemistry.html）

import numpy as np
import scipy.spatial as ssp
from mayavi import mlab
atoms = np.array(
    [[0.0, 0.0, 0.0],
    [0.0, 0.5, 0.5],
    [0.5, 0.0, 0.5],
    [0.5, 0.5, 0.0],
    [0.25, 0.25, 0.25],
    [0.25, 0.75, 0.75],
    [0.75, 0.25, 0.75],
    [0.75, 0.75, 0.25]])
vor = ssp.Voronoi(atoms)

mlab.figure(1, bgcolor=(0, 0, 0), size=(350, 350))
mlab.clf()

# The position of the atoms
atoms_ridge_x = atoms[vor.ridge_points, 0]
atoms_ridge_y = atoms[vor.ridge_points, 1]
atoms_ridge_z = atoms[vor.ridge_points, 2]
tot_uniq = set()
atoms_count = np.zeros(atoms.shape[0], dtype=int)
for i in range(atoms.shape[0]):
    uniq = [j for j, x in enumerate(vor.ridge_points) if i in x]
    uniq = [x for x in uniq if x not in tot_uniq]
    tot_uniq.update(uniq)
    for k in uniq:
        if atoms_count[i] >= 4:
            continue
        mlab.plot3d(atoms_ridge_x[k, :], atoms_ridge_y[k, :], atoms_ridge_z[k, :],
                    tube_radius=0.05, colormap='Reds')
        atoms_count[vor.ridge_points[k]] += 1

H2 = mlab.points3d(atoms[:, 0], atoms[:, 1], atoms[:, 2],
                   scale_factor=0.2,
                   resolution=20,
                   color=(1, 0, 0),
                   scale_mode='none')

mlab.show()

【讨论】：