Making HMR work with dummies

[msuruzhon]

Hello,

I have just discovered that when applying HMR, the dummies are not scaled in accordance with the corresponding hydrogens. I think the expected behaviour would be for the dummies to always have the same mass as the corresponding nondummy atom. How can one make sure this is the case when calling repartitionHydrogenMass()? I am attaching some code for testing purposes.

In addition, it seems that after applying HMR, the dummy gets a negative mass! (-1.00795 g/mol). That would be unexpected behaviour in any case.

Many thanks.

import pytest
import BioSimSpace.Sandpit.Exscientia as BSS

@pytest.fixture
def merged_system():
    ff = "openff_unconstrained-2.0.0"
    mol0 = BSS.Parameters.parameterise("c1ccccc1", ff).getMolecule()
    mol1 = BSS.Parameters.parameterise("c1ccccc1C", ff).getMolecule()
    merged_mol = BSS.Align.merge(mol0, mol1)
    return merged_mol

def test_hmr(merged_system):
    mass0 = merged_system.getAtoms()[12]._sire_object.property("mass0")
    mass1 = merged_system.getAtoms()[12]._sire_object.property("mass1")
    assert mass0 == mass1, "This does not fail"
    merged_system.repartitionHydrogenMass(factor=3)
    BSS.IO.saveMolecules("temp", merged_system, "grotop")
    mass0 = merged_system.getAtoms()[12]._sire_object.property("mass0")
    mass1 = merged_system.getAtoms()[12]._sire_object.property("mass1")
    assert mass0 == mass1, "This will fail, but shouldn't"

[lohedges]

Thanks for this. I hadn't thought about dummies when writing the code and agree with your suggested solution. I'll check what SOMD does too, since I believe I largely followed the implementation there, or at least tested mine against that one when first written.

[lohedges]

I think the issue is that the repartitioning code generates a connectivity object to working the bonding when performing the HMR. This doesn't work properly if an atom is a dummy, since it's based on equilibrium bond radii for specific elements. I think I'll need to adjust this to use the elements from the non-dummy end state instead. Hopefully this should be easy enough to fix.

Thanks for reporting.

[msuruzhon]

Thanks @lohedges, would this work if both endstates have dummies? I think that's the common use case.

[lohedges]

Yes, assuming the velocity conversion is okay I think capping with a warning would be an okay solution.

[lohedges]

I think there are two approaches:

• Include dummies in the repartitioning by using the element property from the opposite end state. In this case dummies are part of the repartitioning so the mass of the end state is consistent before and after repartitioning, but the mass of a dummy atom isn't necessarily the same between end states as it was before.
• Exclude dummies from the repartitioning, then set their mass to that of the opposite end state after repartitioning. In this case, the overall molecule mass of the end state may differ before / after repartitioning (if the sum includes dummies) but the dummy mass is consistent between end states.

[msuruzhon]

Yes I was thinking about that as well, it seems to me that the second approach makes more sense, as one would not expect conservation of mass in an alchemical perturbation setting anyway as you are adding dummies regardless.

[lohedges]

Yes, I agree. I'll check how things are implemented in SOMD and see if @jmichel80 has any input. It's an easy fix to make once we agree on the correct approach. (It's almost like you want to do the partitioning before the merge, i.e. just add the dummies afterwards.)

[lohedges]

SOMD uses the second approach above, i.e. excluding dummy atoms and replacing with the repartitioned mass at the other end state. I'll implement this and check that things are working as expected.

[lohedges]

Okay, this is now fixed, i.e. I get the same behaviour as SOMD, which agrees with the approach agreed upon above. This requires a new build of Sire, so I'll submit a PR when that's done. (Likely Monday.)

Thanks again for catching this.

[msuruzhon]

Many thanks for that @lohedges, I am off until Jan 3rd, so I will only be able to test it then. Happy holidays!

[lohedges]

Thanks, you too. Have a good break!