How to input 3D coordinates from xyz file and connectivity from SMILES in rdkit?

Question

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.

Answer 1

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))

Answer 2

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.