10
$\begingroup$

I am working on a QSAR project where the 3D structural descriptors are an input to a machine learning model. I am generating the descriptors using the python Mordred API (which uses rdkit).

from rdkit import Chem
from mordred import Calculator,descriptors
import pandas as pd

data = pd.read_csv('output_data.csv') # contains SMILES string of all molecules

calc = Calculator(descriptors,ignore_3D=False)

for index,row in data.iterrows():
    mol = Chem.MolFromSmiles(row['SMILES']) # get the SMILES string from each row
    # I need to put in the 3D coordinates from an xyz into 'mol' object here
    descs = list(calc(mol))
    # 
    # ... processing descriptors...
    #

The problem is that if I use rdkit to generate 3D coordinates from the SMILES, it will generate slightly different structures every time. However, I have optimized structure for each of those molecules as xyz files from the calculation that I ran on Orca. I wish to use the optimized structures for the 3D descriptors generation.

How do I do this? The SMILES has the connectivity data, but the xyz has the coordinate data.

I am open to using OpenBabel or any other tool that can be used from Python or the command line, but the main point is to combine connectivity of SMILES with the coordinates of xyz.

$\endgroup$
2
  • $\begingroup$ How is the connectivity stored in an RDKit mol object? Is there anyway to access it once you have made the object? I don't know much about RDKit, but googling around didn't lead to a way to set the connectivity when creating a mol object. $\endgroup$
    – Tyberius
    Dec 10, 2021 at 22:38
  • $\begingroup$ @Tyberius I am not really sure how the connectivity is stored and whether it can be accessed. The documentation of rdkit is not very detailed unfortunately. That is why I asked this question here, in case someone has done this before. I am also open to using other tools apart from rdkit by the way. $\endgroup$
    – S R Maiti
    Dec 10, 2021 at 23:28

2 Answers 2

7
$\begingroup$

Your best bet is to keep the molecule (e.g., as a .mol or .sdf file) around and then simply update the coordinates from the .xyz file as you wish.

As you mention, reading bond orders from XYZ is imperfect - the file simply doesn't contain them and you have to guess. There is no perfect way to do that. Even a highly accurate algorithm for bond order perception .. maybe 97% is still going to make 3% errors.

So rather than attempting to reconstruct the topology, keep it and update the coordinates from the geometry optimization. No loss of information.

You can do this either in Open Babel or RDKit (or other toolkits). Since you mentioned RDKit, here's some rough code that should do that.

m = Chem.MolFromMolFile(sdf_file)

# this assumes whatever program you use doesn't re-order atoms
#  .. which is usually a safe assumption
#  .. so we don't bother tracking atoms
atomic_symbols = []
xyz_coordinates = []

with open(xyz_file, "r") as file:
    for line_number,line in enumerate(file):
        if line_number == 0:
            num_atoms = int(line)
        elif line_number == 1:
            comment = line # might have useful information
        else:
            atomic_symbol, x, y, z = line.split()
            atomic_symbols.append(atomic_symbol)
            xyz_coordinates.append([float(x),float(y),float(z)])

# from https://github.com/rdkit/rdkit/issues/2413
from rdkit.Geometry import Point3D
conf = m.GetConformer()

# in principal, you should check that the atoms match
for i in range(m.GetNumAtoms()):
   x,y,z = xyz_coordinates[i]
   conf.SetAtomPosition(i,Point3D(x,y,z))
$\endgroup$
3
  • $\begingroup$ Thanks! Do you know if rdkit keeps the order of atoms intact when reading SMILES into mol object and then converting it to a .mol file? Otherwise the atoms won't match as you wrote in the answer. $\endgroup$
    – S R Maiti
    Dec 22, 2021 at 10:43
  • 1
    $\begingroup$ In our case, we use RDKit ETDKG to create the initial geometry (and the initial 3D .mol file) so yes, the atoms will match. Well, some versions of the CREST program from the Grimme group do atom reordering, which was annoying. If we're using Open Babel later, we use --gen3d to create the initial .mol file and again, it will match atom order. $\endgroup$ Dec 22, 2021 at 15:31
  • 1
    $\begingroup$ The key thing is "why throw away data" - keep the topology (in whatever format) and just update the coordinates. $\endgroup$ Dec 22, 2021 at 15:31
4
$\begingroup$

The representation of a molecular structure as SMILES string is a (very) reduced one, often with only implicit hydrogen atoms. Recreating a 3D molecular structure with either OpenBabel or RDKit not only adds the atoms together, it will perform a quick/computational affordable force field optimization with some seed, too. (The seed may be randomized, e.g., RDKit manual.) Thus the variation in the atomic coordinates you observe. OpenBabel and RDKit are not written to substitute optimizing an organic structure with high-level theory (like DFT). They may suggest a different conformer than the one your computations with ORCA suggest/experimental findings from docking experiments and structure elucidation (X-ray diffraction, NMR, etc.).

The .xyz file only stores information about atom coordinates. If you want to retain information both about the spatial arrangement as well as connectivity of the atoms in a molecule, a .mol file is better (or, if you want to store multiple models in one container, each about one molecule, a .sdf) with a dedicated connectivity table (see., e.g., Wikipedia's entry, and the current documentation if you want/need to go into detail). RDKit is capable to work with them, too (see e.g., entry from RDKit's documentation), both with the elder v2000 as well as the more recent v3000 format (referring to RDKit version 2021.09.2) in reading and writing.

Check if ORCA can provide the optimized geometry as .mol, or .sdf for your construction of a QSAR.

$\endgroup$
2
  • $\begingroup$ Unfortunately, Orca cannot handle .mol or .sdf files, it can only handle .xyz files. And when OpenBabel converts .xyz to .mol or .sdf, the bond orders it infers are often wrong, which just won't do for fingerprinting. And as I am doing 3D fingerprints, it is better to have the DFT optimized 3D coordinates than the auto generated ones. That is why this is a problem. $\endgroup$
    – S R Maiti
    Dec 12, 2021 at 12:15
  • $\begingroup$ In 2020, Jan Jensen presented xyz2mol as an alternative to OpenBabel (for organic/organometallic molecules) at a RDKit user group meeting (video recording, the public GitHub repository). I'm unable to quantify the differences for the .mol generation by the two programs (OpenBabel and xyz2mol), yet it may be worth to test the later (e.g., in an instance of Miniconda) since it was written for future use in RDKit in mind. $\endgroup$
    – Buttonwood
    Dec 12, 2021 at 13:05

You must log in to answer this question.

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