Skip to content

Small atom mapping debugging #740

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
+42 −26
Open

Small atom mapping debugging #740

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9b8ddb2
Allow non-resonance species to pass isomorphism check
kfir4444 Apr 21, 2024
f37cbf3
Added option to map hydrogen atoms on linear torsions
kfir4444 Apr 21, 2024
24dc585
Minor: torchani -> ob
kfir4444 May 6, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion arc/job/adapters/scripts/ob_script.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@ def save_output_file(path,

def main():
"""
Run a job with torchani guesses.
Run a job with openbabel guesses.
"""
args = parse_command_line_arguments()
try:
Expand Down
66 changes: 41 additions & 25 deletions arc/mapping/engine.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
from arc.species import ARCSpecies
from arc.species.conformers import determine_chirality
from arc.species.converter import compare_confs, sort_xyz_using_indices, translate_xyz, xyz_from_data, xyz_to_str
from arc.species.vectors import calculate_angle, calculate_dihedral_angle, calculate_distance, get_delta_angle
from arc.species.vectors import calculate_angle, calculate_dihedral_angle, calculate_distance, get_delta_angle, VectorsError
from numpy import unique

if TYPE_CHECKING:
Expand Down Expand Up @@ -856,17 +856,23 @@ def map_hydrogens(spc_1: ARCSpecies,
heavy_atom_1_index = atoms_1.index(heavy_atom_1)
for rotor_dict in spc_1.rotors_dict.values():
if heavy_atom_1_index in [rotor_dict['torsion'][1], rotor_dict['torsion'][2]]:
atom_map = add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(
spc_1=spc_1,
spc_2=spc_2,
heavy_atom_1=heavy_atom_1,
heavy_atom_2=heavy_atom_2,
torsion=rotor_dict['torsion'],
atom_map=atom_map,
find_torsion_end_to_map=True,
)
success = True
break
try:
atom_map = add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(
spc_1=spc_1,
spc_2=spc_2,
heavy_atom_1=heavy_atom_1,
heavy_atom_2=heavy_atom_2,
torsion=rotor_dict['torsion'],
atom_map=atom_map,
find_torsion_end_to_map=True,
)
success = True
break
except VectorsError:
# Trying to map a linear torsion
# but the torsion is undefined, since v1 x v2 = 0 if v1 || v2
# using other internal coordinates to map the hydrogens.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

does this fix solve the issue or skips it? We could consider using a different atom instead of the last atom in the torsion

continue
# 3. Check for a pseudo-torsion (may involve multiple bonds) with heavy_atom_1 as a pivot.
if not success:
pseudo_torsion = list()
Expand All @@ -883,17 +889,23 @@ def map_hydrogens(spc_1: ARCSpecies,
pseudo_torsion = [atoms_1.index(atom) for atom in [hydrogen_1, heavy_atom_1, atom_1_3, atom_1_4]]
break
if len(pseudo_torsion):
atom_map = add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(
spc_1=spc_1,
spc_2=spc_2,
heavy_atom_1=heavy_atom_1,
heavy_atom_2=heavy_atom_2,
torsion=pseudo_torsion[::-1],
atom_map=atom_map,
find_torsion_end_to_map=False,
)
success = True
break
try:
atom_map = add_adjacent_hydrogen_atoms_to_map_based_on_a_specific_torsion(
spc_1=spc_1,
spc_2=spc_2,
heavy_atom_1=heavy_atom_1,
heavy_atom_2=heavy_atom_2,
torsion=pseudo_torsion[::-1],
atom_map=atom_map,
find_torsion_end_to_map=False,
)
success = True
break
except VectorsError:
# Trying to map a linear torsion
# but the torsion is undefined, since v1 x v2 = 0 if v1 || v2
# using other internal coordinates to map the hydrogens.
continue
# 4. Check by angles and bond lengths (search for 2 consecutive heavy atoms).
if not success:
atom_1_3, angle_1, bond_length_1 = None, None, None
Expand Down Expand Up @@ -1117,8 +1129,12 @@ def r_cut_p_cut_isomorphic(reactant, product):
bool: ``True`` if they are isomorphic, ``False`` otherwise.
"""
res1 = generate_resonance_structures(object_=reactant.mol.copy(deep=True), save_order = True)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would recommend to modify the generate_resonance_structures func so it always returns at least the original Molecule. Please see where else this func is used and whether this change will be appropriate

for res in res1:
if res.fingerprint == product.mol.fingerprint or product.mol.is_isomorphic(res, save_order=True):
if res1 is not None:
for res in res1:
if res.fingerprint == product.mol.fingerprint or product.mol.is_isomorphic(res, save_order=True):
return True
else:
if reactant.mol.fingerprint == product.mol.fingerprint or product.mol.is_isomorphic(reactant.mol.copy(deep=True), save_order=True):
return True
return False

Expand Down