Skip to content

Commit

Permalink
fix another logical problem: regular Hs (isotope type = 0)
Browse files Browse the repository at this point in the history
not considered for isotope swap
  • Loading branch information
rwxayheee committed Jan 10, 2025
1 parent ad3a81a commit 0906daf
Showing 1 changed file with 14 additions and 18 deletions.
32 changes: 14 additions & 18 deletions meeko/utils/rdkitutils.py
Original file line number Diff line number Diff line change
Expand Up @@ -167,9 +167,6 @@ def get_h_isotope_data(mol: Chem.Mol) -> dict:
parent_atom = copy_mol.GetAtomWithIdx(parent_idx)
# get avail Hs of parent atom
kids_all = [nei for nei in parent_atom.GetNeighbors() if nei.GetAtomicNum() == 1]
# get just enough regular Hs
num_kids_needed = len(isotope_data['parent_idx'][parent_idx]['kids'])
kids_all = kids_all[:num_kids_needed]
# get regular Hs idx
kids_idxs_all = [nei.GetIdx() for nei in kids_all]

Expand All @@ -191,10 +188,7 @@ def get_h_isotope_data(mol: Chem.Mol) -> dict:
# get a list of isotope types on current parent
kids_types = [nei.GetIsotope() for nei in kids_all]
# in case the chirality can't be inverted due to unsolvable problems
if any(len(lst) <= 1 for lst in (
isotope_data['parent_idx'][parent_idx]['kids'], # not enough isotope
set(kids_types), # not enough isotope type
)):
if len(set(kids_types)) <= 1: # not enough isotopes
raise RuntimeError(
f"Unable to recover original chirality by manipulating H sotope positions: \n"
f"Atom # ({parent_idx}) \n"
Expand All @@ -204,14 +198,16 @@ def get_h_isotope_data(mol: Chem.Mol) -> dict:
"or become ambiguous to Chem.rdmolops.AssignStereochemistry due to strained geometry. "
)
# swap assignment of max and min isotope, re-evaluate CIPCode
max_idx = kids_types.index(max(kids_types))
min_idx = kids_types.index(min(kids_types))
copy_mol.GetAtomWithIdx(kids_idxs_all[max_idx]).SetIsotope(kids_types[min_idx])
copy_mol.GetAtomWithIdx(kids_idxs_all[min_idx]).SetIsotope(kids_types[max_idx])
# find the max and min isotope type values
max_isotope_type, min_isotope_type = (max(kids_types), min(kids_types))
# get kiw (kid's index with) max and min isotope types
kiw_max_isotope = kids_idxs_all[kids_types.index(max(kids_types))]
kiw_min_isotope = kids_idxs_all[kids_types.index(min(kids_types))]
copy_mol.GetAtomWithIdx(kiw_max_isotope).SetIsotope(min_isotope_type)
copy_mol.GetAtomWithIdx(kiw_min_isotope).SetIsotope(max_isotope_type)
# re-evaluate current CIPCode from 3D
Chem.AssignStereochemistryFrom3D(copy_mol)
current_cip_code = parent_atom.GetProp('_CIPCode') if parent_atom.HasProp('_CIPCode') else None
current_cip_code = parent_atom.GetProp('_CIPCode') if parent_atom.HasProp('_CIPCode') else None
# a single swap of max and min isotopes normally inverts the chirality, in case not
if current_cip_code!=expected_cip_code:
raise RuntimeError(
Expand All @@ -223,12 +219,12 @@ def get_h_isotope_data(mol: Chem.Mol) -> dict:
"or become ambiguous to Chem.rdmolops.AssignStereochemistry due to strained geometry. "
)
# apply the swap
max_pos = copy_conf.GetAtomPosition(kids_idxs_all[max_idx])
min_pos = copy_conf.GetAtomPosition(kids_idxs_all[min_idx])
H_isotope_idx_with_maxpos = isotope_data['parent_idx'][parent_idx]['kids'][max_idx]
assigned_isotope_pos[H_isotope_idx_with_maxpos] = min_pos
H_isotope_idx_with_minpos = isotope_data['parent_idx'][parent_idx]['kids'][min_idx]
assigned_isotope_pos[H_isotope_idx_with_minpos] = max_pos
max_pos = copy_conf.GetAtomPosition(kiw_max_isotope)
min_pos = copy_conf.GetAtomPosition(kiw_min_isotope)
assigned_isotope_pos[kiw_max_isotope] = min_pos
# only assign positions if min isotope type is greater than 0 (not regular Hs)
if min_isotope_type>0:
assigned_isotope_pos[kiw_min_isotope] = max_pos

return assigned_isotope_pos

Expand Down

0 comments on commit 0906daf

Please sign in to comment.