如何使用pymatgen获得一个结构中的所有键角？

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