Pymatgen structure to rdkit molecule

Question

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?

Answer 1

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.