将组合的索引与值相关联
我正在做一个程序，我需要组合原子之间的距离，或者3D空间中的各个点。下面是一个例子：
文件“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。我这样做是为了计算Pauling电荷，这里你需要看一个原子，确定它有多少键(原子小于2.2安培)，然后看看附着在那个原子上的原子，看看有多少个原子附着在这些原子上。这是非常令人沮丧的，但这一切都将取决于对每个原子的跟踪，而不仅仅是它们的组合。数组可能非常有用。
回答 1
Stack Overflow用户
发布于 2016-06-23 22:05:50
在我们开始之前，让我注意，如果是晶体(我有点怀疑你不是在处理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是可以接受的(所以我选择后者)。
对近邻的后处理使我们脱离了“矮胖”的范畴(你应该尽可能地坚持“矮胖”，这通常是最快的)。对于每个原子，你都想得到一个靠近它的原子的列表。嗯，这不是一个每个原子都有恒定长度的对象，所以你不能把它很好地存储在一个数组中。逻辑的结论是使用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中找到距离足够小的索引，并且一维索引数组存储在生成的list neighbslist的kth元素中。然后检查每个原子的这些数组的长度就很简单了，给出了您的“neihbours数”列表。请注意，我们可以将np.where的输出转换为列表中的list，这样就可以完全保留numpy，但是在下一行中我们必须使用len(neighbs)而不是neighbs.size。
因此，您有两个关键变量，准确地说是两个列表；nearnum[k]是原子k的“近”邻居数( k在range(allpos.shape[0])中，neighbslist[k]是一个一维数字数组，列出了原子k的近邻索引)，所以neighbslist[k][j] ( j in range(nearnum[k]))是range(allpos.shape[0])中的一个数字，不等于k。想想看，这个数组列表结构可能有点难看，所以您可能应该在构建过程中将这个对象转换为一个适当的列表列表(即使这意味着一些开销)。
最后，我只注意到输入数据在文件中。不用担心，那也可以很容易的读到用蒙皮！假设这些空行不在输入名称test中，则可以调用
allpos = np.loadtxt('test',usecols=(1,2,3))将位置矩阵读入变量中。usecols选项允许numpy忽略数据的第一列，该列不是数字列，将导致问题。反正我们也不需要这个。
https://stackoverflow.com/questions/37936527
复制相似问题
腾讯云开发者
