4
$\begingroup$

Context

I am trying to convert a pymatgen bonded structure, such as one generated by CrystalNN to an rdkit molecule.

So I obtain the pymatgen bonded structure with CrystalNN() as shown below:

from pymatgen.io.cif import CifParser
from pymatgen.analysis.local_env import CrystalNN

parser = CifParser('/path/to/cif.cif')
structure = parser.get_structures()[0]

nn = CrystalNN()
adj_matrix = nn.get_bonded_structure(structure)

I can then loop over all the sites and bonds and add them using RWMol from rdkit. I am able to do that as shown below:

from rdkit import Chem

mol = Chem.RWMol() 

#Add atoms
for idx, site in enumerate(structure):
    element = site.as_dict()['species'][0]['element']
    coord = site.coords
    
    atom = Chem.Atom(element)  
    atom.SetDoubleProp('x', coord[0])  
    atom.SetDoubleProp('y', coord[1])  
    atom.SetDoubleProp('z', coord[2])  

    mol.AddAtom(atom) 


#Add bonds
nninfo = nn.get_all_nn_info(structure)

for idx, nbrs in enumerate(nninfo):
    for nidx, nbr in enumerate(nbrs):
        
        nbridx = nninfo[idx][nidx]['site_index']

        try:
            mol.AddBond(idx, nbridx)    
        except:
            continue;

Questions

But this leaves two issues unresolved:

  1. How to limit this just to the molecule, not the entire structure?
  2. How to find the bond-order from CrystalNN?
$\endgroup$
1
  • 1
    $\begingroup$ You have already solved question 2 in your code, you are setting the atomic coordinates using: atom.SetDoubleProp('x', coord[0]) , similarly for y and z. So could you please edit your question accordingly if you are looking for something else? $\endgroup$ May 5 at 16:42

1 Answer 1

3
$\begingroup$

Since you have not mentioned the CIF ID of your crystal, I took a random crystal structure from Crystallography Open Database with CIF ID: 2008748 with which I will answer both your questions 1 and 2.

Initalise the crystal structure:

from pymatgen.io.cif import CifParser
from pymatgen.analysis.local_env import CrystalNN

parser = CifParser('/content/drive/MyDrive/input_files/CIF/2008748.cif')
structure = parser.get_structures()[0]

nn = CrystalNN()
adj_matrix = nn.get_bonded_structure(structure)

Answer for Question 1:

The get_nn_info method returns a list of dictionaries, with each dictionary containing information about a nearest neighbor. The weight key of each dictionary contains the bond order of that neighbor, which you can extract and print.

Code:

nninfo = nn.get_all_nn_info(structure)

for idx, nbrs in enumerate(nninfo):
    for nidx, nbr in enumerate(nbrs):
        bond_order = nninfo[idx][nidx]['weight']
        print(f"The bond order is {bond_order}" )

Answer for Question 2:

This crystal structure has a water molecule. So to extract information related to only the water molecule and not the full structure, you can do as follows:

from pymatgen.core.structure import Molecule
# Define the water molecule
water = Molecule(["O", "H", "H"], [[0.0000, 0.0000, 0.0000], [0.7572, 0.5864, 0.0000], [-0.7572, 0.5864, 0.0000]])

# Find the water molecule in the structure
for i, site in enumerate(structure):
    if site.species_string == "O" and all(site.distance_from_point(point) < 10 for point in water.cart_coords):
        
        water_center = i
        break

# Extract the water molecule from the structure
water_indices = [water_center] + [nn.index for nn in structure.get_neighbors(structure[water_center], 1.5)]
print(water_indices)

You can then use just the molecule for your further analysis.

Note: I obtained the value 10 by initially checking the distance of the water molecule from the site and then plugged it in.

$\endgroup$
2
  • $\begingroup$ Thank you. My question regarding the moelcule was slightly different. How can I get all the different types of molecules in a structure without knowing which molecules might be present in advance? $\endgroup$
    – R Walser
    May 7 at 8:00
  • $\begingroup$ Also, I think your statement that the 'weight' is the bond order is incorrect. The algorithm can return fractional weights too if you set the weighted_cn flag. $\endgroup$
    – R Walser
    May 15 at 6:21

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .