2

所以，使用 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

所以我有一个邻居数组，我也知道它们对应的是哪个原子，现在我只需要获取所有相邻原子的角度。

任何帮助将不胜感激。

4

1 回答 1

0

更新：

所以，我想出了一种方法来做到这一点。我实际上将每个相邻原子转换为它们的笛卡尔坐标。然后，我减去相邻原子所连接的中心原子的坐标。最后，我通过计算每两个向量的点积除以它们的大小的乘积来找到角度，然后取这些值的反余弦值。我使用的代码如下。这可能不是最优雅的做事方式，但它确实完成了工作。如果有人有改进，请随时发表评论。

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
于 2017-02-15T02:35:50.500 回答