Get unique bonds / bond lengths

Is there a way to get just the set of unique bonds / bond lengths? Additionally for a bonded structure, is there an algorithm to output the associated set of bond angles?

I’m trying to achieve a set that looks like

site1, site2, distance

and

site1, site2, site3, angle

or even just the list of indices forming the unique set of bonds which could be mapped e.g. to Structure.get_distance. This information nominally appears in the pymatgen.analysis.graphs.StructureGraph class but I’m struggling to get from this object to the desired output.

Thanks,
Peter

Is there a way to get just the set of unique bonds / bond lengths?

You’re correct that StructureGraph has this functionality in pymatgen.analysis.graphs. This works in two steps,

(1) you create a graph using StructureGraph.with_local_env_strategy() and a suitable strategy from pymatgen.analysis.local_env – these strategies are what determine whether two atoms are bonded, CrystalNN is the usual strategy we recommend.

and (2) call the properties types_and_weights_of_connections, types_of_coordination_environments and weight_statistics might be helpful.

Additionally for a bonded structure, is there an algorithm to output the associated set of bond angles?

There have been a few implementations of this, let me see if I can find a recommended way for you to get this information.

The types_and_weights_of_connections isn’t quite what I was after. I’m trying to use pymatgen as a preprocessing environment for generating input in a structure refinement code. What I’d like to spit out is a list like

# i, j, k, angle
[(2, 0, 10, 117.76368522699826),
 (2, 0, 8, 116.2076335864373),
 (8, 0, 2, 116.20763358643731),
...

so that I can output

# site1 site2 site3 angle
C2 C0 O10 117.76368522699826
C2 C0 O8 116.2076335864373
O8 C0 C2 116.20763358643731
O8 C0 O10 126.00012082009057
O10 C0 C2 117.76368522699826
O10 C0 O8 126.00012082009057
C3 C1 O11 117.76368522699835
C3 C1 O9 116.20763358643727
...

to be parsed as a set of bond angle restraints. Here I’ve added some properties to the Structure instance to generate the unique atom labels.

Below is a simple python approach if there isn’t a more pymatgen-centric approach in existence.

def get_connections(structure):
    """ all connectivity per CrystalNN strategy """
    analyzer = CrystalNN()
    allnn = analyzer.get_all_nn_info(structure)
    rv = []
    for idx, info in enumerate(allnn):
        for el in info:
            rv.append( (idx, el['site_index']) )
    
    return sorted(rv, key=lambda l: l[0])


def get_bond_lengths(structure):
    """ unique bond lengths (i,j) == (j,i) """
    l = []
    for el in get_connections(structure):
        if (not el[::-1] in l) and (not el in l):
            l.append(el)
    rv = []
    for i, j in l:   
        rv.append((i, j, structure.get_distance(i, j)))
    
    return sorted(rv, key=lambda l: l[-1])


def get_bond_angles(structure):
    """ unique bond angles (i,j,k) == (k,j,i) """
    l = []
    con = np.array(get_connections(structure))
    for el in con:
        mask = np.where(con[:, 0] == el[-1])
        for le in con[mask]:
            if not all(el == le[::-1]):
                ijk = (el[0], el[1], le[-1])
                if (not ijk in l) and (not ijk[::-1] in l) :
                    l.append(ijk)
    return sorted([(i, j, k, structure.get_angle(i, j, k) ) for i, j, k in l], key=lambda el: el[1])

Hi @petmetz,

I think this kind of functionality would be useful inside pymatgen. If you wanted, you could open a Pull Request on our GitHub to add this functionality and we can iterate on how it should look there.

You may also want to check out Robocrystallographer, this is a tool built on top of pymatgen and its “intermediate format” contains bond angles I think.

Best,

Matt