Broken dummy atoms logic in SOMD

[jmichel80]

The move to 2023 seems to have broken SOMD. The example below attempts to carry out an RBFE calculation on a protein-ligand system.

I am using feature_pmefep_2023 ; Sire version: 2023.0.0 (28ff974/GitHub/2022-08-25T07:52:15+01:00)

To reproduce the runtime error extract the attached inputs, with that Sire version in your environement

debug.zip

unzip debug.zip 
cd debug/output/lambda-0.00
bash RUNME

Generates on my test computer

(...)
### Running Single Topology Molecular Dynamics Free Energy (v0.2) on node01 ###
###================Setting up calculation=====================###
New run. Loading input and creating restart
lambda is 0.0
Create the System...
Selecting dummy groups
Traceback (most recent call last):
  File "/export/users/julien/miniforge3/envs/sire-dev/share/Sire/scripts/somd-freenrg.py", line 161, in <module>
    OpenMMMD.runFreeNrg(params)
  File "/export/users/julien/miniforge3/envs/sire-dev/lib/python3.9/site-packages/sire/legacy/Tools/__init__.py", line 175, in inner
    retval = func()
  File "/export/users/julien/miniforge3/envs/sire-dev/lib/python3.9/site-packages/sire/legacy/Tools/OpenMMMD.py", line 2275, in runFreeNrg
    system = createSystemFreeEnergy(molecules)
  File "/export/users/julien/miniforge3/envs/sire-dev/lib/python3.9/site-packages/sire/legacy/Tools/OpenMMMD.py", line 1273, in createSystemFreeEnergy
    solute_ref_todummy = solute_ref_todummy.add(
UserWarning: SireError::incompatible_error: The molecules "", number 0, and "LIG", number 1004, are different, and therefore incompatible. (call sire.error.get_last_error_details() for more info)

Manually editing the pert file in debug/input/system.pert to reset all atom types chante to 'du' to a non dummy type resolves the runtime error.

The runtime error is caused by the following line of code

       for x in range(0, ndummies):
            dummy_index = dummies[x].index()
            solute_ref_hard = solute_ref_hard.subtract(
                solute.select(dummy_index)
            )
            solute_ref_todummy = solute_ref_todummy.add(
                solute.select(dummy_index)
            )

@annamherz - could you confirm the provide inputs run without error on a pre-2023 build ?

[jmichel80]

@annamherz note that to run the code in another branch you may have to comment out the line

cutoff type = PME

in the sim.cfg file

I have reproduced the error using the standard cutoffperiodic option so I do not believe the issue is caused by PME specific code changes in that feature branch.

[lohedges]

Thanks, @jmichel80. I should be able to take a look later today.

[lohedges]

Here's a minimal example of the issue using BioSimSpace:

Sire < 2023:

import BioSimSpace as BSS
from Sire.Mol import AtomIdx, Selector_Atom_

# Load an example system and get the first molecule.
mol = BSS.IO.readMolecules("amber/ala/ala*")[0]

# Create an empty selector by selecting all atoms and inverting, as is done by SOMD.
selector0 = m.selectAllAtoms().invert()
print(selector0)
Selector<SireMol::Atom>( [  ] )

# What molecule does this correspond to?
print(selector0.molecule())
Molecule( 2 version 21 : nAtoms() = 22, nResidues() = 3 )

# Create an empty atom selector by constructing the object directly.
selector1 = Selector_Atom()
print(selector1)
Selector<SireMol::Atom>( [  ] )

# What molecule does this correspond to?
print(selector1.molecule())
Molecule( 0 version 0 : nAtoms() = 0, nResidues() = 0 )

# Are these the same?
selector0 == selector1
False

# Add the first molecule to the selector. This works.
selector0.add(mol.getAtoms()[0]._sire_object))
Selector<SireMol::Atom>( [ AtomIdx(0) ] )

# Now try adding to the second selector. This fails with the error you are seeing.
selector1.add(mol.getAtoms()[0]._sire_object))
...
Exception 'SireError::incompatible_error' thrown by the thread 'master'.
The molecules "", number 0, and "", number , are different, and therefore incompatible.

Sire 2023:

import BioSimSpace as BSS
from sire.legacy.Mol import AtomIdx, Selector_Atom_

# Load an example system and get the first molecule.
mol = BSS.IO.readMolecules("amber/ala/ala*")[0]

# Create an empty selector by selecting all atoms and inverting, as is done by SOMD.
selector0 = mol._sire_object.selectAllAtoms().invert()
print(selector0)
Selector<SireMol::Atom>::empty

# What molecule does this correspond to?
print(selector0.molecule())
Molecule( :0      num_atoms=0 num_residues=0 )

# Create an empty atom selector by constructing the object directly.
selector1 = Selector_Atom()
print(selector1)
Selector<SireMol::Atom>::empty

# What molecule does this correspond to?
print(selector1.molecule())
Molecule( :0      num_atoms=0 num_residues=0 )

# Are these the same?
selector0 == selector1
True

# Add the first molecule to the selector. This now fails.
selector0.add(mol.getAtoms()[0]._sire_object))
UserWarning: SireError::incompatible_error: The molecules "", number 0, and "", number 633, are different, and therefore incompatible. (call sire.error.get_last_error_details() for more info)

# Now try adding to the second selector. This still fails.
selector1.add(mol.getAtoms()[0]._sire_object))
...
UserWarning: SireError::incompatible_error: The molecules "", number 0, and "", number 633, are different, and therefore incompatible. (call sire.error.get_last_error_details() for more info)

It looks like an empty selection on a molecule is being converted to a null Selector_Atom_, so is losing reference to the molecule to which the selection was applied.

[jmichel80]

Something doesn't add up since this particular code in SOMD has been around for several years and used for many thousands of RBFEs without this error showing up previously, it's very odd to see a failure on Sire < 2023

[lohedges]

Sorry, it works on < 2023. Apologies if that wasn't clear in the example, or if I made a typo. In my first example, it works if the selector is made from a molecule, as is done in SOMD. For 2023, it fails regardless of how the selector was created.

[lohedges]

Okay, I've made some modifications to SOMD and it now appears to be working. (It starts up and runs without crashing, at least.) Could you check this at your end (perhaps there are other surprises) and if everything is okay then I'll push a fix.

diff --git a/wrapper/Tools/OpenMMMD.py b/wrapper/Tools/OpenMMMD.py
index 74fa1ba3..35c0a0e7 100644
--- a/wrapper/Tools/OpenMMMD.py
+++ b/wrapper/Tools/OpenMMMD.py
@@ -307,10 +307,10 @@ def getSolute(system):
     # matching the perturbed_resnum.val.
 
     # Create the query string.
-    query = f"mol with resnum {perturbed_resnum.val}"
+    query = f"resnum {perturbed_resnum.val}"
 
     # Perform the search.
-    search = system.search(query)
+    search = system.search(query).molecules()
 
     # Make sure there is only one result.
     if len(search) != 1:
@@ -858,6 +858,8 @@ def getDummies(molecule):
 
     from_dummies = None
     to_dummies = None
+    from_non_dummies = []
+    to_non_dummies = []
 
     for x in range(0, natoms):
         atom = atoms[x]
@@ -866,13 +868,17 @@ def getDummies(molecule):
                 from_dummies = molecule.selectAll(atom.index())
             else:
                 from_dummies += molecule.selectAll(atom.index())
-        elif atom.property("final_ambertype") == "du":
+        else:
+            from_non_dummies.append(atom.index())
+        if atom.property("final_ambertype") == "du":
             if to_dummies is None:
                 to_dummies = molecule.selectAll(atom.index())
             else:
                 to_dummies += molecule.selectAll(atom.index())
+        else:
+            to_non_dummies.append(atom.index())
 
-    return to_dummies, from_dummies
+    return to_dummies, from_dummies, to_non_dummies, from_non_dummies
 
 
 def createSystemFreeEnergy(molecules):
@@ -946,10 +952,10 @@ def createSystemFreeEnergy(molecules):
     solute_grp_ref_fromdummy = MoleculeGroup("solute_ref_fromdummy")
 
     solute_ref_hard = solute.selectAllAtoms()
-    solute_ref_todummy = solute_ref_hard.invert()
-    solute_ref_fromdummy = solute_ref_hard.invert()
+    solute_ref_todummy = solute.selectAllAtoms()
+    solute_ref_fromdummy = solute.selectAllAtoms()
 
-    to_dummies, from_dummies = getDummies(solute)
+    to_dummies, from_dummies, to_non_dummies, from_non_dummies = getDummies(solute)
 
     if to_dummies is not None:
         ndummies = to_dummies.count()
@@ -958,7 +964,9 @@ def createSystemFreeEnergy(molecules):
         for x in range(0, ndummies):
             dummy_index = dummies[x].index()
             solute_ref_hard = solute_ref_hard.subtract(solute.select(dummy_index))
-            solute_ref_todummy = solute_ref_todummy.add(solute.select(dummy_index))
+
+    for non_dummy in to_non_dummies:
+        solute_ref_todummy = solute_ref_todummy.subtract(solute.select(non_dummy))
 
     if from_dummies is not None:
         ndummies = from_dummies.count()
@@ -967,7 +975,9 @@ def createSystemFreeEnergy(molecules):
         for x in range(0, ndummies):
             dummy_index = dummies[x].index()
             solute_ref_hard = solute_ref_hard.subtract(solute.select(dummy_index))
-            solute_ref_fromdummy = solute_ref_fromdummy.add(solute.select(dummy_index))
+
+    for non_dummy in from_non_dummies:
+        solute_ref_fromdummy = solute_ref_fromdummy.subtract(solute.select(non_dummy))
 
     solute_grp_ref_hard.add(solute_ref_hard)
     solute_grp_ref_todummy.add(solute_ref_todummy)

The trick was to ensure that we never used an empty selection, which lost reference to the original solute. Instead, we use a full atom selection on the solute and subtract non-dummy atoms in the initial and final state to generate, e.g., solute_ref_todummy. I also needed to update to search string used to extract the perturbable molecule, since this has also changed in 2023. I'll add a single lambda leg FEP test into BioSimSpace so we catch issues like this in future. (SOMD isn't tested as part of Sire itself.)

Cheers.

[lohedges]

I've now added a basic SOMD FEP test to BioSimSpace on my current feature branch here. This fails with the current Sire 2023, but passes with the application of the diff. I'll update the test/s if we find other parts of SOMD that no longer work as expected.

(This does nothing to check that the output is valid, only that it's possible to configure and run a SOMD process without error.)

[lohedges]

I also discovered an issue with the generateDistanceRestraintsDict function which was (incorrectly) assuming that the solute would be MolNum(1). I've switched this to using the getSolute function, which will be correct regardless of the MolNum associated with the solute. The generateDistanceRestraintsDict function was also broken if distance restraints were applied to a single solute system, i.e, it would assume that there were multiple solutes and fail. I've now added error handling for this.

[jmichel80]

I think this feature was not updated in Hannes's branch as work started on a version where the perturbed solute was. Looping @fjclark for awareness as this may become relevant when we migrate the current ABFE work to Sire >=2023.*