[BUG] Asymmetric exclusions for ghost-ghost interactions
Closed this issue · 4 comments
The sire.morph.create_from_pertfile function assumes that the intrascale matrix is the same at lambda = 0 and lambda = 1 and simpy uses the lambda = 0 value for both end states. This leads to issues computing exceptions, leading to an incorrect set of exclusions for the ghost force fields. (Only interactions between ghosts at lambda = 0 are present, those at lambda = 1 are excluded.) BioSimSpace has code to do this, which could possibly be reproduced. However, I think it might not be possible to do this entirely from the lambda = 0 topology and pert file, i.e. some info might be lost.
Extracting the end states and re-merging (with BioSimSpace) seems to recover the correct intrascale matrices. This might be a route that we could take, although it would require extraction of the mapping too, which should be easy enough.
This actually doesn't seem to be limited to create_from_pertfile. Using input that has been merged in both directions, i.e. A-->B and B>>A, I see a different set of exclusions in the ghost-ghost force field, hence the atoms that interact are different depending on the direction of the merge. This isn't resolved by swap-end-states. I'll see if I can cook up a simple example to use for debugging purposes.
Hmm, my debugging script used the incorrect XML file for one merge. I do see this locally for another system, but it's a proprietary molecule. Will keep searching for something that I can post.
Closing since I can't seem to reproduce the error with my original test system. I think this must have been fixed by #237, but I was still testing against one of the old (pre-fix) XML files. I'll re-open if I see it again. The following script shows that (for the given perturbation) things are reversible, both using a regular merged molecule, and one created from a somd1 pert file for the same perturbation:
import BioSimSpace as BSS
import sire as sr
from openmm import XmlSerializer
# Read some ligands from the BioSimSpae tutorial suite.
url = BSS.tutorialUrl()
m0 = BSS.IO.readMolecules([f"{url}/ligand31.prm7.bz2", f"{url}/ligand31.rst7.bz2"])[0]
m1 = BSS.IO.readMolecules([f"{url}/ligand04.prm7.bz2", f"{url}/ligand04.rst7.bz2"])[0]
# Merge with m0 as the reference state.
merged_01 = BSS.Align.merge(m0, m1)
# Merge with m1 as the reference state.
merged_10 = BSS.Align.merge(m1, m0)
# Create two Sire molecules linked to the reference.
mol_01 = sr.morph.link_to_reference(merged_01._sire_object)
mol_10 = sr.morph.link_to_reference(merged_10._sire_object)
map = {
"ghosts_are_light": True,
"check_for_h_by_max_mass": True,
"check_for_h_by_mass": False,
"check_for_h_by_element": False,
"check_for_h_by_ambertype": False,
"fix_perturbable_zero_sigmas": True,
}
d_01 = mol_01.dynamics(perturbable_constraint="h_bonds_not_heavy_perturbed", map=map)
d_10 = mol_10.dynamics(perturbable_constraint="h_bonds_not_heavy_perturbed", map=map)
# Serlialize the OpenMM systems.
with open("debug_01.xml", "w") as f:
f.write(XmlSerializer.serialize(d_01._d._omm_mols.getSystem()))
with open("debug_10.xml", "w") as f:
f.write(XmlSerializer.serialize(d_10._d._omm_mols.getSystem()))
# Now evaulate the atom pairs involved in the ghost-ghost nonbonded force.
pairs_01, _, _ = sr.morph.evaluate_xml_force(mol_01, "debug_01.xml", "ghost-ghost")
pairs_10, _, _ = sr.morph.evaluate_xml_force(mol_10, "debug_10.xml", "ghost-ghost")
print("Pairs 01:")
print(f"{len(pairs_01)} pairs")
for pair in pairs_01:
print(
pair,
pair[0].property("ambertype0"),
pair[1].property("ambertype0"),
pair[0].property("ambertype1"),
pair[1].property("ambertype1"),
)
print("\nPairs 10:")
print(f"{len(pairs_10)} pairs")
for pair in pairs_10:
print(
pair,
pair[0].property("ambertype0"),
pair[1].property("ambertype0"),
pair[0].property("ambertype1"),
pair[1].property("ambertype1"),
)
# Now create SOMD input using the merged molecules.
process01 = BSS.Process.Somd(
merged_01.toSystem(), BSS.Protocol.FreeEnergy(), work_dir="somd_01"
)
process10 = BSS.Process.Somd(
merged_10.toSystem(), BSS.Protocol.FreeEnergy(), work_dir="somd_10"
)
# Recreate the SOMD input using the files.
mols01 = sr.load("somd_01/*7")
mols10 = sr.load("somd_10/*7")
mol_01_somd = sr.morph.create_from_pertfile(mols01[0], "somd_01/somd.pert")
mol_10_somd = sr.morph.create_from_pertfile(mols10[0], "somd_10/somd.pert")
# Create the dynamics objects.
d_01_somd = mol_01_somd.dynamics(
perturbable_constraint="h_bonds_not_heavy_perturbed", map=map
)
d_10_somd = mol_10_somd.dynamics(
perturbable_constraint="h_bonds_not_heavy_perturbed", map=map
)
# Serlialize the OpenMM systems.
with open("debug_01_somd.xml", "w") as f:
f.write(XmlSerializer.serialize(d_01_somd._d._omm_mols.getSystem()))
with open("debug_10_somd.xml", "w") as f:
f.write(XmlSerializer.serialize(d_10_somd._d._omm_mols.getSystem()))
# Now evaulate the atom pairs involved in the ghost-ghost nonbonded force.
pairs_01_somd, _, _ = sr.morph.evaluate_xml_force(
mol_01_somd, "debug_01_somd.xml", "ghost-ghost"
)
pairs_10_somd, _, _ = sr.morph.evaluate_xml_force(
mol_10_somd, "debug_10_somd.xml", "ghost-ghost"
)
print("\nPairs 01 SOMD:")
print(f"{len(pairs_01_somd)} pairs")
for pair in pairs_01_somd:
print(
pair,
pair[0].property("ambertype0"),
pair[1].property("ambertype0"),
pair[0].property("ambertype1"),
pair[1].property("ambertype1"),
)
print("\nPairs 10 SOMD:")
print(f"{len(pairs_10_somd)} pairs")
for pair in pairs_10_somd:
print(
pair,
pair[0].property("ambertype0"),
pair[1].property("ambertype0"),
pair[0].property("ambertype1"),
pair[1].property("ambertype1"),
)This gives:
Pairs 01:
4 pairs
(Atom( H39:42 [ 14.80, 23.18, 18.34] ), Atom( H43:46 [ 17.01, 21.06, 21.15] )) du du hc hc
(Atom( H39:42 [ 14.80, 23.18, 18.34] ), Atom( H44:47 [ 17.41, 22.81, 21.22] )) du du hc hc
(Atom( H40:43 [ 14.63, 21.43, 18.54] ), Atom( H43:46 [ 17.01, 21.06, 21.15] )) du du hc hc
(Atom( H40:43 [ 14.63, 21.43, 18.54] ), Atom( H44:47 [ 17.41, 22.81, 21.22] )) du du hc hc
Pairs 10:
4 pairs
(Atom( H39:39 [ 14.77, 24.13, 19.31] ), Atom( H43:43 [ 17.09, 22.39, 22.29] )) hc hc du du
(Atom( H39:39 [ 14.77, 24.13, 19.31] ), Atom( H44:44 [ 17.66, 24.06, 21.95] )) hc hc du du
(Atom( H40:40 [ 14.47, 22.49, 19.89] ), Atom( H43:43 [ 17.09, 22.39, 22.29] )) hc hc du du
(Atom( H40:40 [ 14.47, 22.49, 19.89] ), Atom( H44:44 [ 17.66, 24.06, 21.95] )) hc hc du du
Pairs 01 SOMD:
4 pairs
(Atom( DU04:42 [ 14.80, 23.18, 18.34] ), Atom( DU08:46 [ 17.01, 21.06, 21.15] )) du du hc hc
(Atom( DU04:42 [ 14.80, 23.18, 18.34] ), Atom( DU09:47 [ 17.41, 22.81, 21.22] )) du du hc hc
(Atom( DU05:43 [ 14.63, 21.43, 18.54] ), Atom( DU08:46 [ 17.01, 21.06, 21.15] )) du du hc hc
(Atom( DU05:43 [ 14.63, 21.43, 18.54] ), Atom( DU09:47 [ 17.41, 22.81, 21.22] )) du du hc hc
Pairs 10 SOMD:
4 pairs
(Atom( H39:39 [ 14.77, 24.13, 19.31] ), Atom( H43:43 [ 17.09, 22.39, 22.29] )) hc hc du du
(Atom( H39:39 [ 14.77, 24.13, 19.31] ), Atom( H44:44 [ 17.66, 24.06, 21.95] )) hc hc du du
(Atom( H40:40 [ 14.47, 22.49, 19.89] ), Atom( H43:43 [ 17.09, 22.39, 22.29] )) hc hc du du
(Atom( H40:40 [ 14.47, 22.49, 19.89] ), Atom( H44:44 [ 17.66, 24.06, 21.95] )) hc hc du du