将组合的索引与值相关联

Question

我正在开发一个程序，我需要组合原子之间的距离或 3D 空间中的各个点。这是一个例子：

文件“测试”包含以下信息：

Ti 1.0 1.0 1.0

O 0.0 2.0 0.0

O 0.0 0.0 0.0

Ti 1.0 3.0 4.0

O 2.0 5.0 0.0

我希望我的代码能够计算点之间距离的所有组合（我已经完成了！），然后，我需要计算一个原子与另一个原子之间的距离小于 2.2 的次数。

这用词令人困惑，所以我将向您展示我目前所掌握的内容。

#!/usr/bin/env python
import sys, math, scipy, itertools
import numpy as np

try:
    infile = sys.argv[1]

except:
    print "Needs file name"
    sys.exit(1)

#opening files for first part
ifile = open(infile, 'r')
coordslist = []

#Creating a file of just coordinates that can be 'mathed on'
for line in ifile:
    pair = line.split()
    atom = (pair[0]); x = float(pair[1]); y = float(pair[2]); z = float(pair[3])
    coordslist += [(x,y,z)]
ifile.close()

#Define distance
def distance(p0,p1):
    return math.sqrt((p0[0] - p1[0])**2 + (p0[1] - p1[1])**2 + (p0[2] - p1[2])**                                          2)

#Initializing for next section
dislist = []
bondslist = []

#Compute distances between all points 1-2, 1-3, 1-4, etc.
for p0, p1 in itertools.combinations(coordslist,2):
    print p0, p1, distance(p0,p1)
    dislist += [distance(p0, p1)]
    if distance(p0,p1) < 2.2:
        bondslist += [(p0, distance(p0,p1))]
print bondslist
print dislist

我不确定列出这些列表是否对我有帮助。到目前为止，他们还没有。

输出是：

(1.0, 1.0, 1.0) (0.0, 2.0, 0.0) 1.73205080757

(1.0, 1.0, 1.0) (0.0, 0.0, 0.0) 1.73205080757

(1.0, 1.0, 1.0) (1.0, 3.0, 4.0) 3.60555127546

(1.0, 1.0, 1.0) (2.0, 5.0, 0.0) 4.24264068712

(0.0, 2.0, 0.0) (0.0, 0.0, 0.0) 2.0

(0.0, 2.0, 0.0) (1.0, 3.0, 4.0) 4.24264068712

(0.0, 2.0, 0.0) (2.0, 5.0, 0.0) 3.60555127546

(0.0, 0.0, 0.0) (1.0, 3.0, 4.0) 5.09901951359

(0.0, 0.0, 0.0) (2.0, 5.0, 0.0) 5.38516480713

(1.0, 3.0, 4.0) (2.0, 5.0, 0.0) 4.58257569496

[((1.0, 1.0, 1.0), 1.7320508075688772), ((1.0, 1.0, 1.0), 1.7320508075688772), ((0.0, 2.0, 0.0), 2.0)]

[1.7320508075688772, 1.7320508075688772, 3.605551275463989, 4.242640687119285, 2.0, 4.242640687119285, 3.605551275463989, 5.0990195135927845, 5.385164807134504, 4.58257569495584]

我需要从这个输出中得到每个原子距离小于 2.2 的次数，例如：

1 2 (because atom 1 has two distances less than 2.2 associated with it)

2 2

3 2 

4 0

5 0

我也需要看看是什么两个原子使小于 2.2 的距离。我这样做是为了计算鲍林费用；这是您需要查看一个原子的地方，确定它有多少个键（距离小于 2.2 埃的原子），然后查看 附着 到该原子的原子，看看有多少原子附在那些上。这非常令人沮丧，但这一切都将取决于跟踪每个原子，而不仅仅是它们的组合。数组可能会非常有用。

我已经检查了here 和here 的帮助，我认为我需要以某种方式结合这些方法。非常感谢任何帮助！

Answer 1

在我们开始之前，让我注意，如果是晶体（我有点怀疑你不是在处理 Ti2O3 分子），你应该小心周期性边界条件，即离所有人较远的最后两个原子可能更接近相邻单元格中的原子。

如果您知道要使用哪些工具，那么您尝试做的事情就非常简单。您正在寻找一种方法来告诉您集合中所有点之间的成对距离。准确地说，执行此操作的函数称为pdist、scipy.spatial.distance.pdist。这可以计算任意维度的任意点集的成对距离，具有任意距离。在您的具体情况下，默认的欧几里得距离就可以了。

一组点的成对矩阵距离（元素[i,j] 告诉您点i 和j 之间的距离）在构造上是对称的，对角线为零。由于这个原因，pdist 的通常实现只返回对角线一侧的非对角线元素，scipy 的版本也不例外。但是，有一个方便的scipy.spatial.distance.squareform 函数可以将包含这种压缩版本的纯非对角对称矩阵的数组转换为满。从那里可以很容易地进行后期处理。

这是我要做的：

import numpy as np
import scipy.spatial as ssp

# atoms and positions:
# Ti 1.0 1.0 1.0
# O  0.0 2.0 0.0
# O  0.0 0.0 0.0
# Ti 1.0 3.0 4.0
# O  2.0 5.0 0.0

# define positions as m*n array, where n is the dimensionality (3)
allpos = np.array([[1.,1,1],  # 1. is lazy for dtype=float64
                   [0,2,0], 
                   [0,0,0],
                   [1,3,4],
                   [2,5,0]])

# compute pairwise distances
alldist_condensed = ssp.distance.pdist(allpos)       # vector of off-diagonal elements on one side
alldist = ssp.distance.squareform(alldist_condensed) # full symmetric distance matrix

# set diagonals to nan (or inf) to avoid tainting our output later
fancy_index = np.arange(alldist.shape[0])
alldist[fancy_index,fancy_index] = np.nan

# find index of "near" neighbours
thresh = 2.2
neighbslist = [np.where(alldist[k,:]<thresh)[0] for k in range(alldist.shape[0])]  # the k'th element is an array containing the indices which are "close" to atom number k

# find total number of "near" neighbours
nearnum = [neighbs.size for neighbs in neighbslist] # the k'th element is the number of atoms which are "close" to atom number k

因此，对于您的具体情况，alldist 包含完整的距离矩阵：

array([[        nan,  1.73205081,  1.73205081,  3.60555128,  4.24264069],
       [ 1.73205081,         nan,  2.        ,  4.24264069,  3.60555128],
       [ 1.73205081,  2.        ,         nan,  5.09901951,  5.38516481],
       [ 3.60555128,  4.24264069,  5.09901951,         nan,  4.58257569],
       [ 4.24264069,  3.60555128,  5.38516481,  4.58257569,         nan]])

如您所见，我手动将对角线元素设置为np.nan。这是必要的，因为我打算检查该矩阵中小于thresh 的元素，并且对角线中的零肯定符合条件。在我们的例子中，np.inf 对于这些元素来说同样是一个不错的选择，但是如果您想要获得比thresh 彼此之间更远 的点怎么办？显然对于这种情况-np.inf 或np.nan 是可以接受的（所以我选择了后者）。

近邻的后处理使我们脱离了 numpy 的领域（你应该尽可能地坚持使用 numpy，这通常是最快的）。对于每个原子，您想要获取靠近它的那些原子的列表。好吧，这不是每个原子都具有恒定长度的对象，因此您不能将其很好地存储在数组中。合乎逻辑的结论是使用list，但是你可以使用所有python并使用列表推导来构造这个列表（上面的提醒）：

neighbslist = [np.where(alldist[k,:]<thresh)[0] for k in range(alldist.shape[0])]  # the k'th element is an array containing the indices which are "close" to atom number k

这里np.where 将在k 行中找到距离足够小的索引，并且一维索引数组存储在结果列表neighbslist 的第kth 元素中。然后检查每个原子的这些数组的长度是微不足道的，为您提供“附近的邻居数”列表。请注意，我们可以将 np.where 的输出转换为列表 comp 中的 list 以完全保留 numpy，但是我们将不得不在下一行中使用 len(neighbs) 而不是 neighbs.size。

所以，你有两个关键变量，准确地说是两个列表； nearnum[k] 是原子 k 的“近”邻居数（k 在 range(allpos.shape[0]) 中，neighbslist[k] 是一个一维 numpy 数组，列出了原子 k 的近索引，所以 neighbslist[k][j]（对于j in range(nearnum[k])) 是range(allpos.shape[0]) 中的一个数字，不等于k。想想看，这个数组列表结构可能有点难看，所以你应该把这个对象转换为在构建过程中正确的列表列表（即使这意味着一些开销）。

---

我最后才注意到您的输入数据在文件中。不用担心，也可以使用 numpy 轻松读取！假设那些空行不在你输入的名字test中，你可以调用

allpos = np.loadtxt('test',usecols=(1,2,3))

将位置矩阵读入变量。 usecols 选项让numpy 忽略数据的第一列，这不是数字，会导致问题。反正我们真的不需要那个。