Skip to content

Fix inconsistencies in merge code #488

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

Merged
lohedges merged 2 commits into from
Dec 5, 2025
Merged

Fix inconsistencies in merge code #488

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
13d0fc0
Fix inconsistencies in merge code.
lohedges Dec 5, 2025
2d5043c
Make sure not to connect zeroed bonds.
lohedges Dec 5, 2025
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
238 changes: 115 additions & 123 deletions python/BioSimSpace/Align/_merge.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ def merge(
merged : Sire.Mol.Molecule
The merged molecule.
"""
from sire.legacy import CAS as _SireCAS
from sire.legacy import Mol as _SireMol
from sire.legacy import Base as _SireBase
from sire.legacy import MM as _SireMM
Expand Down Expand Up @@ -1027,45 +1028,29 @@ def merge(
"or 'allow_ring_size_change' options."
)

# Create the connectivity object.
conn = _SireMol.Connectivity(edit_mol.info()).edit()

# Connectivity in the merged molecule.
conn0 = _SireMol.Connectivity(edit_mol.info()).edit()
conn1 = _SireMol.Connectivity(edit_mol.info()).edit()

bonds_atoms0 = {
frozenset([bond.atom0(), bond.atom1()])
for bond in edit_mol.property("bond0").potentials()
}
bonds_atoms1 = {
frozenset([bond.atom0(), bond.atom1()])
for bond in edit_mol.property("bond1").potentials()
}

for atom0, atom1 in bonds_atoms0:
conn0.connect(atom0, atom1)

for atom0, atom1 in bonds_atoms1:
conn1.connect(atom0, atom1)

# We only add a bond to the total connectivity if it's defined in both states
# This results in a "broken" topology if one writes it in GROMACS
# But GROMACS can't handle bond breaks anyway
for atom0, atom1 in bonds_atoms0 & bonds_atoms1:
conn.connect(atom0, atom1)

conn = conn.commit()
# Connect the bonded atoms in both end states.
for bond in edit_mol.property("bond0").potentials():
# Only connect bonds with non-zero force constants.
ab = _SireMM.AmberBond(bond.function(), _SireCAS.Symbol("r"))
if ab.k() != 0.0:
conn0.connect(bond.atom0(), bond.atom1())
conn0 = conn0.commit()
for bond in edit_mol.property("bond1").potentials():
# Only connect bonds with non-zero force constants.
ab = _SireMM.AmberBond(bond.function(), _SireCAS.Symbol("r"))
if ab.k() != 0.0:
conn1.connect(bond.atom0(), bond.atom1())
conn1 = conn1.commit()

# Get the connectivity of the two molecules.
# Get the original connectivity of the two molecules.
c0 = molecule0.property("connectivity")
c1 = molecule1.property("connectivity")

# Check that the merge hasn't modified the connectivity.

# The checking was blocked when merging a protein
if roi is None:
# molecule0
for x in range(0, molecule0.num_atoms()):
Expand All @@ -1083,7 +1068,7 @@ def merge(
idy_map = mol0_merged_mapping[idy]

# Was a ring opened/closed?
is_ring_broken = _is_ring_broken(c0, conn, idx, idy, idx_map, idy_map)
is_ring_broken = _is_ring_broken(c0, conn1, idx, idy, idx_map, idy_map)

# A ring was broken and it is not allowed.
if is_ring_broken and not allow_ring_breaking:
Expand All @@ -1095,7 +1080,7 @@ def merge(

# Did a ring change size?
is_ring_size_change = _is_ring_size_changed(
c0, conn, idx, idy, idx_map, idy_map
c0, conn1, idx, idy, idx_map, idy_map
)

# A ring changed size and it is not allowed.
Expand All @@ -1113,7 +1098,7 @@ def merge(
)

# The connectivity has changed.
if c0.connection_type(idx, idy) != conn.connection_type(
if c0.connection_type(idx, idy) != conn1.connection_type(
idx_map, idy_map
):
# The connectivity changed for an unknown reason.
Expand Down Expand Up @@ -1141,7 +1126,7 @@ def merge(
idy_map = mol1_merged_mapping[idy]

# Was a ring opened/closed?
is_ring_broken = _is_ring_broken(c1, conn, idx, idy, idx_map, idy_map)
is_ring_broken = _is_ring_broken(c1, conn0, idx, idy, idx_map, idy_map)

# A ring was broken and it is not allowed.
if is_ring_broken and not allow_ring_breaking:
Expand All @@ -1153,7 +1138,7 @@ def merge(

# Did a ring change size?
is_ring_size_change = _is_ring_size_changed(
c1, conn, idx, idy, idx_map, idy_map
c1, conn0, idx, idy, idx_map, idy_map
)

# A ring changed size and it is not allowed.
Expand All @@ -1171,7 +1156,7 @@ def merge(
)

# The connectivity has changed.
if c1.connection_type(idx, idy) != conn.connection_type(
if c1.connection_type(idx, idy) != conn0.connection_type(
idx_map, idy_map
):
# The connectivity changed for an unknown reason.
Expand Down Expand Up @@ -1266,90 +1251,95 @@ def merge(
iterlen = len(molecule1_roi_idx)
iterrange = molecule1_roi_idx

for x in range(0, iterlen):
for x in range(0, iterlen):
# Convert to an AtomIdx.
idx = iterrange[x]
idx = _SireMol.AtomIdx(idx)

# Map the index to its position in the merged molecule.
idx_map = mol1_merged_mapping[idx]

for y in range(x + 1, iterlen):
idy = iterrange[y]
# Convert to an AtomIdx.
idx = iterrange[x]
idx = _SireMol.AtomIdx(idx)
idy = _SireMol.AtomIdx(idy)

# Map the index to its position in the merged molecule.
idx_map = mol1_merged_mapping[idx]

for y in range(x + 1, iterlen):
idy = iterrange[y]
# Convert to an AtomIdx.
idy = _SireMol.AtomIdx(idy)
idy_map = mol1_merged_mapping[idy]

# Map the index to its position in the merged molecule.
idy_map = mol1_merged_mapping[idy]
conn_type = conn0.connection_type(idx_map, idy_map)
# The atoms aren't bonded.
if conn_type == 0:
clj_scale_factor = _SireMM.CLJScaleFactor(1, 1)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)

conn_type = conn0.connection_type(idx_map, idy_map)
# The atoms aren't bonded.
if conn_type == 0:
clj_scale_factor = _SireMM.CLJScaleFactor(1, 1)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)
# The atoms are part of a dihedral.
elif conn_type == 4:
clj_scale_factor = _SireMM.CLJScaleFactor(
ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor()
)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)

# The atoms are part of a dihedral.
elif conn_type == 4:
clj_scale_factor = _SireMM.CLJScaleFactor(
ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor()
)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)
# The atoms are bonded
else:
clj_scale_factor = _SireMM.CLJScaleFactor(0, 0)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)

# The atoms are bonded
else:
clj_scale_factor = _SireMM.CLJScaleFactor(0, 0)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)
# Now copy in all intrascale values from molecule0 into clj_nb_pairs0.

# Now copy in all intrascale values from molecule0 into both
# clj_nb_pairs matrices.
# Work out the iteration length and range.
if roi is None:
iterlen = molecule0.num_atoms()
iterrange = list(range(molecule0.num_atoms()))
# When region of interest is defined, perfrom loop from these indices
else:
for roi_res in roi:
# Get atom indices of ROI residue in molecule0
# Get atom indices of ROI residue in molecule0.
molecule0_roi = molecule0.residues()[roi_res]
molecule0_roi_idx = [a.index() for a in molecule0_roi]
iterlen = len(molecule0_roi_idx)
iterrange = molecule0_roi_idx

# Perform a triangular loop over atoms from molecule0.
for x in range(0, iterlen):
# Convert to an AtomIdx.
idx = iterrange[x]
idx = _SireMol.AtomIdx(idx)
# Perform a triangular loop over atoms from molecule0.
for x in range(0, iterlen):
# Convert to an AtomIdx.
idx = iterrange[x]
idx = _SireMol.AtomIdx(idx)

# Map the index to its position in the merged molecule.
idx_map = mol0_merged_mapping[idx]

for y in range(x + 1, iterlen):
idy = iterrange[y]
# Convert to an AtomIdx.
idy = _SireMol.AtomIdx(idy)

# Map the index to its position in the merged molecule.
idy_map = mol0_merged_mapping[idy]

conn_type = conn0.connection_type(idx_map, idy_map)
# The atoms aren't bonded.
if conn_type == 0:
clj_scale_factor = _SireMM.CLJScaleFactor(1, 1)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)

# The atoms are part of a dihedral.
elif conn_type == 4:
clj_scale_factor = _SireMM.CLJScaleFactor(
ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor()
)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)
# Map the index to its position in the merged molecule.
idx_map = mol0_merged_mapping[idx]

# The atoms are bonded
else:
clj_scale_factor = _SireMM.CLJScaleFactor(0, 0)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)
for y in range(x + 1, iterlen):
idy = iterrange[y]
# Convert to an AtomIdx.
idy = _SireMol.AtomIdx(idy)

# Map the index to its position in the merged molecule.
idy_map = mol0_merged_mapping[idy]

# Get the connection type between the atoms.
conn_type = conn0.connection_type(idx_map, idy_map)

# The atoms aren't bonded.
if conn_type == 0:
clj_scale_factor = _SireMM.CLJScaleFactor(1, 1)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)

# The atoms are part of a dihedral.
elif conn_type == 4:
clj_scale_factor = _SireMM.CLJScaleFactor(
ff.electrostatic14_scale_factor(),
ff.vdw14_scale_factor(),
)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)

# The atoms are bonded
else:
clj_scale_factor = _SireMM.CLJScaleFactor(0, 0)
clj_nb_pairs0.set(idx_map, idy_map, clj_scale_factor)

# Finally, copy the intrascale from molecule1 into clj_nb_pairs1.

# Work out the iteration length and range.
if roi is None:
iterlen = molecule1.num_atoms()
iterrange = list(range(molecule1.num_atoms()))
Expand All @@ -1361,40 +1351,42 @@ def merge(
iterlen = len(molecule1_roi_idx)
iterrange = molecule1_roi_idx

# Perform a triangular loop over atoms from molecule1.
for x in range(0, iterlen):
# Convert to an AtomIdx.
idx = iterrange[x]
idx = _SireMol.AtomIdx(idx)
# Perform a triangular loop over atoms from molecule1.
for x in range(0, iterlen):
# Convert to an AtomIdx.
idx = iterrange[x]
idx = _SireMol.AtomIdx(idx)

# Map the index to its position in the merged molecule.
idx = mol1_merged_mapping[idx]
# Map the index to its position in the merged molecule.
idx = mol1_merged_mapping[idx]

for y in range(x + 1, iterlen):
idy = iterrange[y]
# Convert to an AtomIdx.
idy = _SireMol.AtomIdx(idy)
for y in range(x + 1, iterlen):
idy = iterrange[y]
# Convert to an AtomIdx.
idy = _SireMol.AtomIdx(idy)

# Map the index to its position in the merged molecule.
idy = mol1_merged_mapping[idy]
# Map the index to its position in the merged molecule.
idy = mol1_merged_mapping[idy]

conn_type = conn1.connection_type(idx, idy)
# Get the connection type between the atoms.
conn_type = conn1.connection_type(idx, idy)

if conn_type == 0:
clj_scale_factor = _SireMM.CLJScaleFactor(1, 1)
clj_nb_pairs1.set(idx, idy, clj_scale_factor)
if conn_type == 0:
clj_scale_factor = _SireMM.CLJScaleFactor(1, 1)
clj_nb_pairs1.set(idx, idy, clj_scale_factor)

# The atoms are part of a dihedral.
elif conn_type == 4:
clj_scale_factor = _SireMM.CLJScaleFactor(
ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor()
)
clj_nb_pairs1.set(idx, idy, clj_scale_factor)
# The atoms are part of a dihedral.
elif conn_type == 4:
clj_scale_factor = _SireMM.CLJScaleFactor(
ff.electrostatic14_scale_factor(),
ff.vdw14_scale_factor(),
)
clj_nb_pairs1.set(idx, idy, clj_scale_factor)

# The atoms are bonded
else:
clj_scale_factor = _SireMM.CLJScaleFactor(0, 0)
clj_nb_pairs1.set(idx, idy, clj_scale_factor)
# The atoms are bonded
else:
clj_scale_factor = _SireMM.CLJScaleFactor(0, 0)
clj_nb_pairs1.set(idx, idy, clj_scale_factor)

# Store the two molecular components.
edit_mol.set_property("molecule0", molecule0)
Expand Down
Loading