Skip to content

Commit 14e376c

Browse files
lohedgeslohedges
authored
Merge pull request #313 from OpenBioSim/fix_312
Fix issue #312
2 parents ede60ba + edadeee commit 14e376c

File tree

5 files changed

+374
-225
lines changed

5 files changed

+374
-225
lines changed

‎python/BioSimSpace/Align/_align.py

Lines changed: 161 additions & 111 deletions
Original file line numberDiff line numberDiff line change
@@ -1547,22 +1547,47 @@ def _rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_ma
15471547
mol0 = molecule0._getSireObject()
15481548
mol1 = molecule1._getSireObject()
15491549

1550+
# Do we have a monatomic ion?
1551+
if (molecule0.nAtoms() == 1) or (molecule1.nAtoms() == 1):
1552+
is_ion = True
1553+
else:
1554+
is_ion = False
1555+
15501556
# Convert the mapping to AtomIdx key:value pairs.
15511557
sire_mapping = _to_sire_mapping(mapping)
15521558

15531559
# Perform the alignment, mol0 to mol1.
1554-
try:
1560+
if len(mapping) == 1 and is_ion:
1561+
idx0 = list(mapping.keys())[0]
1562+
idx1 = list(mapping.values())[0]
1563+
# Replace the coordinates of the mapped atom with those of the reference.
15551564
mol0 = (
1556-
mol0.move().align(mol1, _SireMol.AtomResultMatcher(sire_mapping)).molecule()
1565+
mol0.edit()
1566+
.atom(idx0)
1567+
.setProperty(
1568+
property_map0.get("coordinates", "coordinates"),
1569+
mol1.atom(idx1).property(
1570+
property_map1.get("coordinates", "coordinates")
1571+
),
1572+
)
1573+
.molecule()
1574+
.commit()
15571575
)
1558-
except Exception as e:
1559-
msg = "Failed to align molecules based on mapping: %r" % mapping
1560-
if "Could not calculate the single value decomposition" in str(e):
1561-
msg += ". Try minimising your molecular coordinates prior to alignment."
1562-
if _isVerbose():
1563-
raise _AlignmentError(msg) from e
1564-
else:
1565-
raise _AlignmentError(msg) from None
1576+
else:
1577+
try:
1578+
mol0 = (
1579+
mol0.move()
1580+
.align(mol1, _SireMol.AtomResultMatcher(sire_mapping))
1581+
.molecule()
1582+
)
1583+
except Exception as e:
1584+
msg = "Failed to align molecules based on mapping: %r" % mapping
1585+
if "Could not calculate the single value decomposition" in str(e):
1586+
msg += ". Try minimising your molecular coordinates prior to alignment."
1587+
if _isVerbose():
1588+
raise _AlignmentError(msg) from e
1589+
else:
1590+
raise _AlignmentError(msg) from None
15661591

15671592
# Return the aligned molecule.
15681593
return _Molecule(mol0)
@@ -2509,6 +2534,12 @@ def _score_rdkit_mappings(
25092534
.commit()
25102535
)
25112536

2537+
# Do we have a monatomic ion?
2538+
if (molecule0.nAtoms() == 1) or (molecule1.nAtoms() == 1):
2539+
is_ion = True
2540+
else:
2541+
is_ion = False
2542+
25122543
# Get the set of matching substructures in each molecule. For some reason
25132544
# setting uniquify to True removes valid matches, in some cases even the
25142545
# best match! As such, we set uniquify to False and account for duplicate
@@ -2575,61 +2606,68 @@ def _score_rdkit_mappings(
25752606
break
25762607

25772608
if is_valid:
2578-
# Rigidly align molecule0 to molecule1 based on the mapping.
2579-
if scoring_function == "RMSDALIGN":
2580-
try:
2581-
molecule0 = (
2582-
molecule0.move()
2583-
.align(
2584-
molecule1, _SireMol.AtomResultMatcher(sire_mapping)
2585-
)
2586-
.molecule()
2587-
)
2588-
except Exception as e:
2589-
if (
2590-
"Could not calculate the single value decomposition"
2591-
in str(e)
2592-
):
2593-
is_gsl_error = True
2594-
gsl_exception = e
2595-
else:
2596-
msg = (
2597-
"Failed to align molecules when scoring based on mapping: %r"
2598-
% mapping
2609+
# If there is only a single atom in the mapping and one molecule
2610+
# has one atom, e.g. an ion, then skip the alignment.
2611+
if len(mapping) == 1 and is_ion:
2612+
mappings.append(mapping)
2613+
scores.append(0.0)
2614+
else:
2615+
# Rigidly align molecule0 to molecule1 based on the mapping.
2616+
if scoring_function == "RMSDALIGN":
2617+
try:
2618+
molecule0 = (
2619+
molecule0.move()
2620+
.align(
2621+
molecule1,
2622+
_SireMol.AtomResultMatcher(sire_mapping),
2623+
)
2624+
.molecule()
25992625
)
2600-
if _isVerbose():
2601-
raise _AlignmentError(msg) from e
2626+
except Exception as e:
2627+
if (
2628+
"Could not calculate the single value decomposition"
2629+
in str(e)
2630+
):
2631+
is_gsl_error = True
2632+
gsl_exception = e
26022633
else:
2603-
raise _AlignmentError(msg) from None
2604-
# Flexibly align molecule0 to molecule1 based on the mapping.
2605-
elif scoring_function == "RMSDFLEXALIGN":
2606-
molecule0 = flexAlign(
2607-
_Molecule(molecule0),
2608-
_Molecule(molecule1),
2609-
mapping,
2610-
property_map0=property_map0,
2611-
property_map1=property_map1,
2612-
)._sire_object
2613-
2614-
# Append the mapping to the list.
2615-
mappings.append(mapping)
2616-
2617-
# We now compute the RMSD between the coordinates of the matched atoms
2618-
# in molecule0 and molecule1.
2619-
2620-
# Initialise lists to hold the coordinates.
2621-
c0 = []
2622-
c1 = []
2623-
2624-
# Loop over each atom index in the map.
2625-
for idx0, idx1 in sire_mapping.items():
2626-
# Append the coordinates of the matched atom in molecule0.
2627-
c0.append(molecule0.atom(idx0).property("coordinates"))
2628-
# Append the coordinates of atom in molecule1 to which it maps.
2629-
c1.append(molecule1.atom(idx1).property("coordinates"))
2630-
2631-
# Compute the RMSD between the two sets of coordinates.
2632-
scores.append(_SireMaths.getRMSD(c0, c1))
2634+
msg = (
2635+
"Failed to align molecules when scoring based on mapping: %r"
2636+
% mapping
2637+
)
2638+
if _isVerbose():
2639+
raise _AlignmentError(msg) from e
2640+
else:
2641+
raise _AlignmentError(msg) from None
2642+
# Flexibly align molecule0 to molecule1 based on the mapping.
2643+
elif scoring_function == "RMSDFLEXALIGN":
2644+
molecule0 = flexAlign(
2645+
_Molecule(molecule0),
2646+
_Molecule(molecule1),
2647+
mapping,
2648+
property_map0=property_map0,
2649+
property_map1=property_map1,
2650+
)._sire_object
2651+
2652+
# Append the mapping to the list.
2653+
mappings.append(mapping)
2654+
2655+
# We now compute the RMSD between the coordinates of the matched atoms
2656+
# in molecule0 and molecule1.
2657+
2658+
# Initialise lists to hold the coordinates.
2659+
c0 = []
2660+
c1 = []
2661+
2662+
# Loop over each atom index in the map.
2663+
for idx0, idx1 in sire_mapping.items():
2664+
# Append the coordinates of the matched atom in molecule0.
2665+
c0.append(molecule0.atom(idx0).property("coordinates"))
2666+
# Append the coordinates of atom in molecule1 to which it maps.
2667+
c1.append(molecule1.atom(idx1).property("coordinates"))
2668+
2669+
# Compute the RMSD between the two sets of coordinates.
2670+
scores.append(_SireMaths.getRMSD(c0, c1))
26332671

26342672
# No mappings were found.
26352673
if len(mappings) == 0:
@@ -2732,6 +2770,12 @@ def _score_sire_mappings(
27322770
.commit()
27332771
)
27342772

2773+
# Do we have a monatomic ion?
2774+
if (molecule0.nAtoms() == 1) or (molecule1.nAtoms() == 1):
2775+
is_ion = True
2776+
else:
2777+
is_ion = False
2778+
27352779
# Initialise a list to hold the mappings.
27362780
mappings = []
27372781

@@ -2751,54 +2795,60 @@ def _score_sire_mappings(
27512795
break
27522796

27532797
if is_valid:
2754-
# Rigidly align molecule0 to molecule1 based on the mapping.
2755-
if scoring_function == "RMSDALIGN":
2756-
try:
2757-
molecule0 = (
2758-
molecule0.move()
2759-
.align(molecule1, _SireMol.AtomResultMatcher(mapping))
2760-
.molecule()
2761-
)
2762-
except Exception as e:
2763-
msg = (
2764-
"Failed to align molecules when scoring based on mapping: %r"
2765-
% mapping
2766-
)
2767-
if _isVerbose():
2768-
raise _AlignmentError(msg) from e
2769-
else:
2770-
raise _AlignmentError(msg) from None
2771-
# Flexibly align molecule0 to molecule1 based on the mapping.
2772-
elif scoring_function == "RMSDFLEXALIGN":
2773-
molecule0 = flexAlign(
2774-
_Molecule(molecule0),
2775-
_Molecule(molecule1),
2776-
_from_sire_mapping(mapping),
2777-
property_map0=property_map0,
2778-
property_map1=property_map1,
2779-
)._sire_object
2780-
2781-
# Append the mapping to the list.
2782-
mapping = _from_sire_mapping(mapping)
2783-
mapping = dict(sorted(mapping.items()))
2784-
mappings.append(mapping)
2785-
2786-
# We now compute the RMSD between the coordinates of the matched atoms
2787-
# in molecule0 and molecule1.
2788-
2789-
# Initialise lists to hold the coordinates.
2790-
c0 = []
2791-
c1 = []
2792-
2793-
# Loop over each atom index in the map.
2794-
for idx0, idx1 in mapping.items():
2795-
# Append the coordinates of the matched atom in molecule0.
2796-
c0.append(molecule0.atom(idx0).property("coordinates"))
2797-
# Append the coordinates of atom in molecule1 to which it maps.
2798-
c1.append(molecule1.atom(idx1).property("coordinates"))
2799-
2800-
# Compute the RMSD between the two sets of coordinates.
2801-
scores.append(_SireMaths.getRMSD(c0, c1))
2798+
# If there is only a single atom in the mapping and one molecule
2799+
# has one atom, e.g. an ion, then skip the alignment.
2800+
if len(mapping) == 1 and is_ion:
2801+
mappings.append(mapping)
2802+
scores.append(0.0)
2803+
else:
2804+
# Rigidly align molecule0 to molecule1 based on the mapping.
2805+
if scoring_function == "RMSDALIGN":
2806+
try:
2807+
molecule0 = (
2808+
molecule0.move()
2809+
.align(molecule1, _SireMol.AtomResultMatcher(mapping))
2810+
.molecule()
2811+
)
2812+
except Exception as e:
2813+
msg = (
2814+
"Failed to align molecules when scoring based on mapping: %r"
2815+
% mapping
2816+
)
2817+
if _isVerbose():
2818+
raise _AlignmentError(msg) from e
2819+
else:
2820+
raise _AlignmentError(msg) from None
2821+
# Flexibly align molecule0 to molecule1 based on the mapping.
2822+
elif scoring_function == "RMSDFLEXALIGN":
2823+
molecule0 = flexAlign(
2824+
_Molecule(molecule0),
2825+
_Molecule(molecule1),
2826+
_from_sire_mapping(mapping),
2827+
property_map0=property_map0,
2828+
property_map1=property_map1,
2829+
)._sire_object
2830+
2831+
# Append the mapping to the list.
2832+
mapping = _from_sire_mapping(mapping)
2833+
mapping = dict(sorted(mapping.items()))
2834+
mappings.append(mapping)
2835+
2836+
# We now compute the RMSD between the coordinates of the matched atoms
2837+
# in molecule0 and molecule1.
2838+
2839+
# Initialise lists to hold the coordinates.
2840+
c0 = []
2841+
c1 = []
2842+
2843+
# Loop over each atom index in the map.
2844+
for idx0, idx1 in mapping.items():
2845+
# Append the coordinates of the matched atom in molecule0.
2846+
c0.append(molecule0.atom(idx0).property("coordinates"))
2847+
# Append the coordinates of atom in molecule1 to which it maps.
2848+
c1.append(molecule1.atom(idx1).property("coordinates"))
2849+
2850+
# Compute the RMSD between the two sets of coordinates.
2851+
scores.append(_SireMaths.getRMSD(c0, c1))
28022852

28032853
# No mappings were found.
28042854
if len(mappings) == 0:

‎python/BioSimSpace/Process/_amber.py

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2856,9 +2856,11 @@ def is_exe(fpath):
28562856
# order their path accordingly, or use the exe keyword argument.
28572857
if results:
28582858
for exe in results:
2859-
exe = _pathlib.Path(p) / exe
2860-
if is_exe(exe):
2861-
return str(exe)
2859+
# Exclude "locally enhanced sampling" executables.
2860+
if "LES" not in exe:
2861+
exe = _pathlib.Path(p) / exe
2862+
if is_exe(exe):
2863+
return str(exe)
28622864

28632865
msg = (
28642866
"'BioSimSpace.Process.Amber' is not supported. "

0 commit comments

Comments
 (0)