Как получить все углы связи внутри конструкции, используя пиматген?

Question

Итак, используя pymatgen, у меня есть объект структуры. То, что я хочу сделать, это получить все углы связи внутри структуры. Я мог бы зациклить каждый атом, чтобы получить все углы связи, но это включало бы каждый атом, независимо от того, насколько далеко они друг от друга, что, конечно, является проблемой.

Теперь я могу найти соседей каждого центрального атома, используя функцию "get_neighbors", однако я не уверен, куда идти, тем более, что функция "get_angle" принимает целочисленные значения, а не объекты атомарного сайта.

Ниже приведен код, который у меня есть:

import pymatgen as mg
import numpy as np


s = mg.Structure.from_file('POSCAR')

atoms = s.atomic_numbers

van = [x for x in atoms if x == 23]
length = len(van)
nb = ['NONE']*length

x = 0

while x < length:

    nb[x] = s.get_neighbors(s[x],2.4)
    x += 1

Итак, у меня есть массив соседей, и я знаю, какому атому они тоже соответствуют, теперь мне просто нужно получить углы для всех соседних атомов.

Любая помощь будет принята с благодарностью.

Answer 1

Обновить:

Итак, я нашел способ сделать это. Я фактически преобразовал каждый соседний атом в их декартовы координаты. Затем я вычел координаты центрального атома, к которому присоединены соседние атомы. Наконец, я нашел углы, рассчитав скалярное произведение каждых двух векторов, деленное на произведение их величин, а затем взяв арккосинус этих значений. Код, который я использовал ниже. Это может быть не самый элегантный способ сделать что-то, но он выполняет свою работу. Если у кого-то есть улучшения, не стесняйтесь комментировать.

import random
import itertools as iter
import pymatgen as mg
import numpy as np
import math


def all_angles(POSCAR,amin,amax):

    s = mg.Structure.from_file(POSCAR)

    atoms = s.atomic_numbers

    van = [x for x in atoms if x == 23]
    length = len(van)
    nb = ['NONE']*length

    x = 0
    n = 0

    while x < length:

        nb[x] = s.get_neighbors(s[x],2.4)
        x += 1

    w = 100
    h = len(van)
    oxygen = [[0 for x in range(w)] for y in range(h)]

    x = 0
    y = 0

    while x < len(nb):

        y = 0
        n = 0

        while y < len(nb[x]):

            oxygen[x][n] = (nb[x][y][0]).coords
            n += 1
            y += 1

        x += 1

    x = 0

    while x < len(oxygen):
        oxygen[x] = [n for n in oxygen[x] if not isinstance(n,int)]
        x += 1

    van = [x for x in range(0,len(van))]

    x = 0
    while x < len(van):
        van[x] = s[van[x]].coords
        x += 1

    van = [np.array(x) for x in van]

    x = 0
    while x < len(van):

        oxygen[x] = [np.subtract(oxygen[x][y],van[x]) for y in range(0,len(oxygen[x]))]
        x += 1

    combo = [[0 for x in range(0,1000)] for y in range(0,len(van))]

    r = 0

    while r < len(van):
        x = 0
        for subset in iter.combinations(oxygen[r],2):
            combo[r][x] = subset
            x += 1
        r += 1

    x = 0
    while x < len(combo):
        combo[x] = [c for c in combo[x] if c != 0]
        x += 1


    angles =  [[0 for x in range(0,1000)] for y in range(0,len(van))]

    x = 0

    while x < len(combo):

        group = combo[x]

        y = 0

        while y < len(group):
            angles[x][y] = math.degrees(math.acos(np.dot(group[y][0],group[y][1])/(np.linalg.norm(group[y][0])*np.linalg.norm(group[y][1]))))
            y += 1

        angles[x] = [round(n,3) for n in angles[x] if n != 0 and n > amin and n < amax]

        x += 1

    angles = np.concatenate(angles,axis=0)    

    return angles