Связывание индексов комбинации со значением
Я работаю над программой, для которой мне нужны комбинации расстояний между атомами или различными точками в трехмерном пространстве. Вот пример:
Файл '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 ангстрем), затем посмотреть на атомы, связанные с этим атомом, и посмотреть, сколько атомов присоединено к ним. Это ужасно расстраивает, но все будет зависеть от отслеживания каждого атома, а не только их комбинации. Массив, вероятно, будет чрезвычайно полезен.
Я проверял здесь и здесь помощь, и я думаю, что мне нужно каким-то образом объединить эти методы. Любая помощь невероятно ценится!
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 для которого расстояние достаточно мало, а массив индексов 1d хранится в k ый элемент результирующего списка neighbslist, Затем тривиально проверять длину этих массивов для каждого атома, давая вам список "число ближайших соседей". Обратите внимание, что мы могли бы привести к выводу np.where к list в списке Comp оставить numy полностью, но тогда нам пришлось бы использовать len(neighbs) вместо neighbs.size в следующей строке.
Итак, у вас есть две ключевые переменные, если быть точным, два списка; nearnum[k] число "ближних" соседей для атома k (с k в range(allpos.shape[0]), а также neighbslist[k] является 1d массивом с перечислением ближайших индексов для атома k, так neighbslist[k][j] (за j в range(nearnum[k])) число в range(allpos.shape[0]) не равно k, Подумайте об этом, эта конструкция списка массивов, вероятно, немного уродлива, поэтому вы, вероятно, должны привести этот объект к правильному списку списков во время построения (даже если это потребует некоторых накладных расходов).
Я только заметил в конце, что ваши входные данные находятся в файле. Не волнуйтесь, это также можно легко прочитать, используя NumPy! Предполагая, что эти пустые строки не в вашем имени входа test, ты можешь позвонить
allpos = np.loadtxt('test',usecols=(1,2,3))
читать матрицу положения в вашу переменную. usecols опция позволяет numpy игнорируйте первый столбец данных, который не является числовым и может вызвать проблемы. Нам все равно это не нужно.