Adjust positions of hydrogen isotopes as if they were regular Hs

[rwxayheee]

This is for #295 and #298.

Added a function set_h_isotope_atom_coords to replace SetTerminalAtomCoords for the assignment of H isotope coordinates with AddHs, and with preserved chirality if from input molecular graph.

More description in comments.

[rwxayheee]

I tried to implement the embed-align method which is a different approach to this PR. The main advantage of doing that is it has an almost seamless way to handle isotopes and isomers. But the re-generated mol conformation can still look different (even the RMSD is low, like 0.5) from the docked pose, which sometimes possesses a strange torsion. With that I don't think it's applicable at the moment. Embedding the isotope-containing fragment is another way, but again if the docked pose is distorted, there's inevitably a shift of the export ligand. Since we also need to recover the exact position of polar hydrogens, we don't want the ligand atom positions to change.

Tomorrow I will go back to the Add-RemoveHs method and continue this (for the chirality rebuild, or check, i'm not sure..). I had many misconceptions about RDKit Smiles handling before this. A lot of what I previously commented on #297 and #295 was wrong. It cost me so much time to realize the way I wrote Smiles has many problems ..

[rwxayheee]

Documenting for future reference:

I tried many different ways to call SetTerminalAtomCoords (on complete mol, on mol with partially assigned conf, etc.) and in different versions. It reproduces AddHs geometry for single atom, when the other atoms' coordinates are reasonable.

However, when multiple atoms (isotopes) are missing a valid coordinates (or just at (0,0,0)), the positions from SetTerminalAtomCoords can go wrong even with an ideal geometry of the heavy atoms. See an example at the end of this comment.

I feel like this is because it's a function intended for single coordinate assignment, assuming the surrounding geometry is complete. This behavior has led to lots of confusion I had these days, but now I think the undesired positions have very little to do with docking poses being distorted. The unwanted geometry results from how the SetTerminalAtomCoords function is designed.

Thus, the real problem might be the usage of SetTerminalAtomCoords and the condition of the molecule conformation when this function is being called. This PR will implement a replacement that's specific for H-isotope (this issue). But in the future, I will keep looking and see if there are other replacements in RDKit.

It should be noted that AddHs does not necessarily return the best geometry, either - there can be close contacts between hydrogens.

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Geometry import Point3D

"""
rdkit.__version__
'2024.09.3'
"""

# SMILES in Issue 295
smi = "[2H]C([2H])([2H])[C@H](O)C([2H])([2H])CCCn1c(=O)c2c(ncn2C)n(C)c1=O"
mol = Chem.MolFromSmiles(smi)
AllChem.EmbedMolecule(mol, AllChem.ETKDG())

mol_with_h = Chem.AddHs(mol, addCoords=True)

# Ideal, as reference
writer = Chem.SDWriter("mol_0.sdf")
writer.write(mol_with_h)
writer.close()

for atom in mol_with_h.GetAtoms():
    idx = atom.GetIdx()
    if atom.GetIsotope()>0:
        neigh = atom.GetNeighbors()
        AllChem.SetTerminalAtomCoords(mol_with_h, idx, neigh[0].GetIdx())

# OK, if isotopes are set one at a time? 
writer = Chem.SDWriter("mol_1.sdf")
writer.write(mol_with_h)
writer.close()

conformer = mol_with_h.GetConformer()
for atom in mol_with_h.GetAtoms():
    idx = atom.GetIdx()
    if atom.GetIsotope()>0:
        conformer.SetAtomPosition(idx, Point3D(0,0,0))

for atom in mol_with_h.GetAtoms():
    idx = atom.GetIdx()
    if atom.GetIsotope()>0:
        neigh = atom.GetNeighbors()
        AllChem.SetTerminalAtomCoords(mol_with_h, idx, neigh[0].GetIdx())

# Not OK - has distorted angels
writer = Chem.SDWriter("mol_2.sdf")
writer.write(mol_with_h)
writer.close()

[diogomart]

I agree: the problem is how I'm using SetTerminalAtomCoords. Based on its source in RDKit's AddHs.cpp, it uses the positions of all atoms in the molecule. We are assigning positions to heavy atoms but leaving deuterium undefined or at (0,0,0). Then, when we call SetTerminalAtomCoords to calculate a position for a given atom, the undefined or zeroed positions of nearby atoms bound to the same heavy atom are taken into account, causing erroneous geometries.

[rwxayheee]

[1a5373b] is a working version. I will do some cleanup, add comments and run more checks tomorrow.

Added a function to replace SetTerminalAtomCoords. It uses AddHs to assign coordinates and tries to correct _CIPCode if the initial assignment does not match with input. This function can be used like a version ofAddHs that includes H isotopes.

Input:
an RDKit mol, it can be under-processing (no regular Hs, new conformer not set yet),
the incomplete conformer,
an optional flag to return fully-processed mol copy (all Hs added, with a complete conformer)

Output:
a dictionary (H isotope atom idx in input mol -> assigned coordinate in Point3D), or a fully-processed mol

Current state:
Can pass existing tests
Output for #295:
DB15122_1a5373b.sdf.txt

The output looks exactly like the output when this PR was at c9435bb.

The major improvement is handling of chirality, and I also refactored and encapsulated the additional process, to insert the minimal amount of codes into rdkit_mol_create.py.

[rwxayheee]

[b1a2835] is a bugfix. I changed the CIP chirality tag assignment methods. Previously, I wasn't using the method that's 3D-aware to re-evaluate the chirality after the initial placement and the isotope swap.

This PR should be ready at this point..! Some comments for future reference:

1. Efficiency and implementation of the added function

The added function, set_h_isotope_atom_coords, is here as a specific replacement to SetTerminalAtomCoords and is encapsulated to be called at the specific location. But because the function uses AddHs, the isotope coordinate assignment could have been deferred to add_hydrogens. I didn't implement the return_mol option, as I don't want not introduce more than necessary changes to the existing RDKitMolCreate procedure. Overall I hope this is ok a temporary fix. When we have a better replacement in RDKit, we can just remove this function.

2. Another occurrence of SetTerminalAtomCoord

There's another occurrence of SetTerminalAtomCoord in molsetup.py. But none of the outer function is currently in use.