[BUG] Asymmetric exclusions for ghost-ghost interactions

[lohedges]

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.

[lohedges]

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.

[lohedges]

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.

[lohedges]

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.

[lohedges]

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