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:
- How to limit this just to the molecule, not the entire structure?
- How to find the bond-order from CrystalNN?
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$