investigate why stereochemistry is not being preserved when embedding conformers with rdkit

[joshkamm]

• [x] See if I can politely ask rdkit to maintain stereocenters - probably yes but nontrivial
• [x] See if I can hardcode an atom as a stereocenter using rdkit - also probably yes but undesirable

I spent some time trying to understand what was going on and looking through the documentation. It seems that rdkit doesn't automatically discover and assign chiral atoms and that this must be done explicitly before calling EmbedMultipleConfs.

It looks like rdkit has recent activity in this area: rdkit/rdkit#5330

A few options going forward:

• See if confab preserves stereochemistry more naturally by rotating torsions from an initial structure and switch to that
  • Check what I did in the C++ version for conformer searching with openbabel
• Use AllChem.AssignStereochemistryFrom3D which appears to set the chiral tags based on rdkit's legacy implementation
  • This doesn't appear to identify the spiro center in Alonso's ligand as chiral
  • I think it would be possible to explicitly label an atom as chiral
• Ask a question in rdkit github discussions about how to use the newer Chem.FindPotentialStereo() implementation to assign chiral centers together with assigning the absolute stereochemistry from a 3d structure

I'm currently trying out openbabel confab and it doesn't appear to have a parallel implementation but at least it seems to naturally preserve stereochemistry. The openbabel python API is built with SWIG which is much less friendly with communicating data structures than rdkit's boost.python.

• [x] Test the updated implementation on a simple example
• [x] Test the updated implementation on a more complex example
• [x] Clean up code and retest
• [x] Look through pull request changes
• [x] Update documentation

[joshkamm]

Test conformer generation code with openbabel and rdkit for future reference:

# EXPERIMENTING WITH RDKIT CONFORMER GENERATION
rdkit_mol = stk_ligand.to_rdkit_mol()
from rdkit import Chem
from rdkit.Chem import AllChem
print(Chem.FindMolChiralCenters(rdkit_mol, includeUnassigned=True, useLegacyImplementation=False))
# AllChem.AssignStereochemistryFrom3D(rdkit_mol, confId=0, replaceExistingTags=True)
si = Chem.FindPotentialStereo(rdkit_mol)
for element in si:
    print(f'  Type: {element.type}, Which: {element.centeredOn}, Specified: {element.specified}, Descriptor: {element.descriptor} ')
pass

# EXPERIMENTING WITH OPEN BABEL CONFORMER GENERATION
from openbabel import pybel as pb
pybel_mol = next(pb.readfile('xyz', str(ligand_path)))
# print(len(pybel_mol.conformers))
cs = pb.ob.OBConformerSearch()
cs.Setup(pybel_mol.OBMol, 10, 5, 5, 5) # numConformers, numChildren, mutability, convergence
cs.Search()
cs.GetConformers(pybel_mol.OBMol)
print(pybel_mol.OBMol.NumConformers())
stk_conformers = []
for i in range(pybel_mol.OBMol.NumConformers()):
    pybel_mol.OBMol.SetConformer(i)
    stk_conformers.append(pybel_mol_to_stk_mol(pybel_mol))
stk_list_to_xyz_file(stk_conformers, 'test_pybel_generation.xyz')    rdkit_mol = stk_ligand.to_rdkit_mol()
    conf_ids = Chem.AllChem.EmbedMultipleConfs(
        rdkit_mol,
        numConfs=config.initial_conformers,
        randomSeed=40,
        pruneRmsThresh=config.initial_rms_threshold,
        numThreads=config.num_cpus
    )
    stk_conformers = [stk_ligand.with_position_matrix(rdkit_mol.GetConformer(conf_id).GetPositions())
                   for conf_id in conf_ids]

Previous code for using rdkit to generate conformers:

    rdkit_mol = stk_ligand.to_rdkit_mol()
    conf_ids = Chem.AllChem.EmbedMultipleConfs(
        rdkit_mol,
        numConfs=config.initial_conformers,
        randomSeed=40,
        pruneRmsThresh=config.initial_rms_threshold,
        numThreads=config.num_cpus
    )
    stk_conformers = [stk_ligand.with_position_matrix(rdkit_mol.GetConformer(conf_id).GetPositions())
                   for conf_id in conf_ids]