How to get all bond angles within a structure, using pymatgen?

Question

So, using pymatgen, I have a structure object. What I want to do is get all of the bond angles within the structure. I could loop over every atom to get all bond angles, but this would include every atom, regardless of how far apart they are, which of course is a problem.

Now, I can find the neighbors of each central atom using the "get_neighbors" function, however, I'm not sure where to go from here, especially because the "get_angle" function takes integer values and not atomic site objects.

Below is the code I have thus far:

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

So I have an array of the neighbors and I know which atom they correspond too, now I just need to get the angles for all neighboring atoms.

Any help would be greatly appreciated.

Answer 1

Update:

So, I've figured out a way of doing this. I actually converted each neighboring atom into their Cartesian coordinates. Then, I subtracted the coordinates of the central atom to which the neighboring atoms are attached. Lastly, I found the angles by calculating the dot product of each two vectors divided by the product of their magnitudes, then taking the arc cosine of these values. The code which I used is below. This may not be the most elegant way of doing things, but it does get the job done. If anyone has improvements, feel free to comment.

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