将组合的索引与值相关联

Question

我正在做一个程序，我需要原子之间的距离的组合，或者三维空间中的不同点。举个例子：

文件“test”包含以下信息：

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，^{}。这可以计算任意维度上任意点集的两两距离，具有任何类型的距离。在您的特定情况下，默认的欧几里德距离就可以了。你知道吗

一组点的成对矩阵距离（元素[i,j]告诉你点i和j之间的距离）通过构造是对称的，对角线上有零。因此，pdist的通常实现只返回对角线一侧的非对角线元素，scipy的版本也不例外。但是，有一个方便的^{}函数，它可以将包含这种纯非对角对称矩阵压缩版本的数组转换为完整的数组。从那里很容易进行后期处理。你知道吗

我会这么做：

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。因为这是必要的，因为这比矩阵的对角线元素要少。在我们的例子中，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行中找到距离足够小的索引，索引的1d数组存储在结果列表neighbslist的第k元素中。然后检查每个原子的这些数组的长度就很简单了，这样就可以得到“附近邻居的数量”列表。请注意，我们可以将np.where的输出强制转换为list comp中的list以完全保留numpy，但这样我们就不得不在下一行中使用len(neighbs)而不是neighbs.size。你知道吗

因此，您有两个关键变量，准确地说是两个列表；nearnum[k]是原子k（其中k位于range(allpos.shape[0])中）的“近邻”数，neighbslist[k]是列出原子k的近邻索引的1d numpy数组，因此neighbslist[k][j]（对于range(nearnum[k])中的j）是range(allpos.shape[0])中不等于k的数。想想看，这个数组列表的构造可能有点难看，所以在构造过程中，您可能应该将这个对象转换为一个适当的列表列表（即使这意味着有点开销）。你知道吗

---

我只注意到最后你的输入数据在一个文件里。不用担心，那也可以用numpy随时阅读！假设这些空行不在您的输入名test中，您可以调用

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

将位置矩阵读入变量。usecols选项允许numpy忽略数据的第一列，这不是数字，并且会导致问题。反正我们也不需要。你知道吗