Сопоставление индексов комбинации со значением

Question

Я работаю над программой, для которой мне нужны комбинации расстояний между атомами или различные точки в трехмерном пространстве. Вот пример:

Файл «тест» содержит следующую информацию: (!, Который я сделал)

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, но тогда вы можете идти питон и использовать список понимание, чтобы построить этот список (напоминание сверху):

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. Тогда тривиально проверять длину этих массивов для каждого атома, предоставляя вам список «число ближайших neihbours». Обратите внимание, что мы могли бы сделать вывод np.where в list в списке comp, чтобы оставить numpy полностью, но тогда нам пришлось бы использовать len(neighbs) вместо neighbs.size в следующей строке.

Итак, у вас есть две ключевые переменные: два списка, чтобы быть точным; nearnum[k] это число «близких» соседей для атома k (с k в range(allpos.shape[0]) и neighbslist[k] является 1d NumPy массив перечисляя рядом индексов для атома k, так neighbslist[k][j] (для j в range(nearnum[k])) является числом в range(allpos.shape[0]) не равна k Подумайте об этом, это построение списков массивов, вероятно, немного уродливое, поэтому вы должны, вероятно, включить этот объект в правильный список списков во время строительства (даже если это означает немного накладных расходов).

---

В конце я заметил, что ваши входные данные находятся в файле. Не волнуйтесь, это также можно легко прочитать, используя NumPy! Если предположить, что эти пустые строки не на ваше имя входного test, вы можете вызвать

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

читать матрицу позиции в переменный. Опция usecols позволяет игнорировать первый столбец данных, который не является числовым, и может вызвать проблемы. В любом случае, нам это действительно не нужно.