Skip to content
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.

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

Small atom mapping debugging #740

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
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
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 @@
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:

Check warning on line 871 in arc/mapping/engine.py

View check run for this annotation

Codecov / codecov/patch

arc/mapping/engine.py#L871

Added line #L871 was not covered by tests
# 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

Check warning on line 875 in arc/mapping/engine.py

View check run for this annotation

Codecov / codecov/patch

arc/mapping/engine.py#L875

Added line #L875 was not covered by tests
# 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 @@
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:

Check warning on line 904 in arc/mapping/engine.py

View check run for this annotation

Codecov / codecov/patch

arc/mapping/engine.py#L904

Added line #L904 was not covered by tests
# 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

Check warning on line 908 in arc/mapping/engine.py

View check run for this annotation

Codecov / codecov/patch

arc/mapping/engine.py#L908

Added line #L908 was not covered by tests
# 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 @@
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
Loading