Tests for Receptor-Ligand Restraints

[fjclark]

Hello,

This accompanies my recent pull request to Sire, and implements a test for Boresch restraints. I've also implemented a test for the multiple distance restraints which are already present in Sire.

I've followed the style of test_combRules.py, copying over most of OpenMMMD.py and using this to set up systems with and without restraints, then testing for an expected energy difference.

Thanks,
Finlay

[lohedges]

Thanks for this. Actually, the approach should be to test the OpenMMD.py installed by Sire, not a duplicate here. That way your test will still be valid if we make further edits to the script, i.e. you aren't just testing a frozen snapshot. I mentioned this for the test_combRules.py PR but it was never updated and I needed to merge things across quickly. (I'll add a note to split this out at some point.) To use functionality from OpenMMD.py you can just import as follows, e.g:

from Sire.Tools.OpenMMMD import *

Just to note that with @chryswoods' feat_web we are beginning the work of bring the tests into the main Sire repository to help keep things in sync. For space reasons, we will also be moving to an approach of keeping test input files elsewhere (on a remote server) and downloading them dynamically during testing. (The new Sire has functionality to do this.) This won't affect this PR. I will just test against your SIre PR to validate that it works and merge when happy. We can then move it into the main repo along with other tests for Sire.Move when the time comes.

Cheers.

[fjclark]

Thanks very much for the comments - I will update this to use functionality from OpenMMMD.py shortly.

Cheers.

[lohedges]

No problem. I'm just trying to update the test_combRules.py file, which could serve as an example.

[lohedges]

Okay, I've re-written the test as follows:

#
# Sire script to compare Sire and SOMD single point energies
#

from nose.tools import assert_almost_equal

from Sire.Tools import OpenMMMD

# Overwrite default parameters.

OpenMMMD.topfile = OpenMMMD.Parameter("topfile", "../io/combRules/toluene_methane.top",
                    """File name of the topology file containing the system to be simulated.""")

OpenMMMD.crdfile = OpenMMMD.Parameter("crdfile", "../io/combRules/toluene_methane.crd",
                    """File name of the coordinate file containing the coordinates of the
                       system to be simulated.""")

OpenMMMD.morphfile = OpenMMMD.Parameter("morphfile", "../io/combRules/toluene_methane.pert",
                      """Name of the morph file containing the perturbation to apply to the system.""")

OpenMMMD.andersen = OpenMMMD.Parameter("thermostat", False,
      """Whether or not to use the Andersen thermostat (needed for NVT or NPT simulation).""")

OpenMMMD.barostat = OpenMMMD.Parameter("barostat", False, """Whether or not to use a barostat (needed for NPT simulation).""")

OpenMMMD.cutoff_type = OpenMMMD.Parameter("cutoff type", "nocutoff", """The cutoff method to use during the simulation.""")

OpenMMMD.constraint = OpenMMMD.Parameter("constraint", "none", """The constraint model to use during dynamics.""")

OpenMMMD.platform = OpenMMMD.Parameter("platform", "Reference", """Which OpenMM platform should be used to perform the dynamics.""")

def test_combining_rules(rule, verbose=False):
    # Set the combining rule.
    OpenMMMD.combining_rules = OpenMMMD.Parameter("combining rules", rule,
                            """Combining rules to use for the non-bonded interactions.""")

    amber = OpenMMMD.Amber()

    (molecules, space) = amber.readCrdTop(OpenMMMD.crdfile.val, OpenMMMD.topfile.val)
    system = OpenMMMD.createSystemFreeEnergy(molecules)

    system = OpenMMMD.setupForceFieldsFreeEnergy(system, space)
    if OpenMMMD.debug_seed.val:
        ranseed = OpenMMMD.debug_seed.val
    else:
        ranseed = OpenMMMD.RanGenerator().randInt(100000, 1000000)

    moves = OpenMMMD.setupMovesFreeEnergy(system, ranseed, 
            OpenMMMD.gpu.val, OpenMMMD.lambda_val.val)

    # Get energy from Sire
    nrg_sire = system.energy()
    # Get energy from SOMD
    mdmoves = moves.moves()[0]
    integrator = mdmoves.integrator()
    nrg_somd = integrator.getPotentialEnergy(system)

    nrg_diff = (nrg_sire - nrg_somd).value()

    if verbose:
        print ("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
        if rule == 'arithmetic':
            print ("TESTING ARITHMETIC COMBINING RULES")
        else:
            print ("TESTING GEOMETRIC COMBINING RULES")
        print ("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")

        print ("* The energy from Sire is: %s" %nrg_sire)
        print ("* The energy from SOMD is: %s" % nrg_somd)

        if rule == 'arithmetic':
            print (" For the arithmetic combining rules the single point energy difference between sire and somd at lambda %s is %s " % (OpenMMMD.lambda_val.val, nrg_diff) )
        else:
            print (" For the geometric combining rules the single point energy difference between sire and somd at lambda %s is %s " % (OpenMMMD.lambda_val.val, nrg_diff) )

    assert_almost_equal(nrg_diff, 0.0, 4)

if __name__ == '__main__':
    test_combining_rules("arithmetic", True)
    test_combining_rules("geometric", True)

Interestingly, the test now fails for geometric combining rules with the following output:

Create the System...
Selecting dummy groups
Creating force fields...
Setting up moves...
Using debugging seed number 756320
Created a MD move that uses OpenMM for all molecules on 0
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
TESTING ARITHMETIC COMBINING RULES
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* The energy from Sire is: 2.68694 kcal mol-1
* The energy from SOMD is: 2.68694 kcal mol-1
 For the arithmetic combining rules the single point energy difference between sire and somd at lambda 0.0 is -6.120802709119744e-06
Create the System...
Selecting dummy groups
Creating force fields...
Setting up moves...
Using debugging seed number 207294
Created a MD move that uses OpenMM for all molecules on 0
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
TESTING GEOMETRIC COMBINING RULES
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* The energy from Sire is: 2.61428 kcal mol-1
* The energy from SOMD is: 2.68694 kcal mol-1
 For the geometric combining rules the single point energy difference between sire and somd at lambda 0.0 is -0.07266638703239758
Traceback (most recent call last):
  File "test_combRules.py", line 80, in <module>
    test_combining_rules("geometric", True)
  File "test_combRules.py", line 76, in test_combining_rules
    assert_almost_equal(nrg_diff, 0.0, 4)
  File "/home/lester/sire.app/lib/python3.7/unittest/case.py", line 906, in assertAlmostEqual
    raise self.failureException(msg)
AssertionError: -0.07266638703239758 != 0.0 within 4 places (0.07266638703239758 difference)

I'll check to see if I've done something wrong. (Perhaps I'm missing a change to one of the default parameters.) If not, this is a good example of why you should test against OpenMMMD.py within Sire, not an extracted and modified version, since updates may have since broken the functionality that you added. (I've re-run the old test and it still passes.)

Looking at the code in the original test vs that in the current OpenMMMD.py it is clear that things are quite different, i.e. for the hydrogen_mass_repartitioning_factor parameter:

test_combRules.py

hydrogen_mass_repartitioning_factor = Parameter("hydrogen mass repartitioning factor",None,
                                     """If not None and is a number, all hydrogen atoms in the molecule will
                                        have their mass increased by the input factor. The atomic mass of the heavy atom
                                        bonded to the hydrogen is decreased to keep the mass constant.""")

OpenMMMD.py

hydrogen_mass_repartitioning_factor = \
    Parameter('hydrogen mass repartitioning factor', 1.0,
              f'If larger than {HMR_MIN} (maximum is {HMR_MAX}), all hydrogen '
              'atoms in the molecule will have their mass increased by this '
              'factor. The atomic mass of the heavy atom bonded to the '
              'hydrogen is decreased to keep the total mass constant '
              '(except when this would lead to a heavy atom to be lighter '
              'than a minimum mass).')

Looking at the git blame for those lines I see:

git blame -L 189,196 wrapper/Tools/OpenMMMD.py 
32bdab798 (Hannes Loeffler 2020-11-18 09:36:20 +0000 189) hydrogen_mass_repartitioning_factor = \
32bdab798 (Hannes Loeffler 2020-11-18 09:36:20 +0000 190)     Parameter('hydrogen mass repartitioning factor', 1.0,
32bdab798 (Hannes Loeffler 2020-11-18 09:36:20 +0000 191)               f'If larger than {HMR_MIN} (maximum is {HMR_MAX}), all hydrogen '
32bdab798 (Hannes Loeffler 2020-11-18 09:36:20 +0000 192)               'atoms in the molecule will have their mass increased by this '
32bdab798 (Hannes Loeffler 2020-11-18 09:36:20 +0000 193)               'factor. The atomic mass of the heavy atom bonded to the '
32bdab798 (Hannes Loeffler 2020-11-18 09:36:20 +0000 194)               'hydrogen is decreased to keep the total mass constant '
32bdab798 (Hannes Loeffler 2020-11-18 09:36:20 +0000 195)               '(except when this would lead to a heavy atom to be lighter '
32bdab798 (Hannes Loeffler 2020-11-18 09:36:20 +0000 196)               'than a minimum mass).')

Those edits are fairly old. I'm wondering if the combining rules in the current version of OpenMMMD.py actually works correctly, since it clearly isn't being tested against the most recent version. I'll try to see if I can track down any more differences. @jmichel80: Do you have any thoughts?

I've also noted that @chryswoods has removed test_combRules.py because of complications porting to the new API. I wonder if there will be similar issues for your test, or whether my modifications would solve this issue?

[lohedges]

Interestingly, the SOMD energies do agree, but the Sire energies do not. Perhaps something needs to be set for the Sire system. I'll look in the functions within the old test script to see if anything has been edited when setting up the system.

[jmichel80]

Hum that's concerning. I can see how something could have slipped by because we are not using in current projects forcefields with geometric combining rules. This could have an impact on QUBE work done currently by @djcole56 and @bieniekmateusz . I would expect the Sire single point energies to change for the same molecule when changing between arithmetic and geometric combining rules so the test likely fails because the test doesn't actually update SOMD to use geometric combining rules for the second set of single point energy calculations

[djcole56]

Thanks for the heads up, but I believe we've switched all our FF work to arithmetic combining rules now so shouldn't affect any current projects.

[fjclark]

Thanks for the example @lohedges - I've updated my script to match.

[lohedges]

Great. I just overrode the same parameters that were edited in the test_combRules.py script, although some had now effect. You'll likely need to tweak whichever are relevant to your own tests, so delete any that aren't needed.

Cheers.

[fjclark]

Thanks. I've already made sure the overridden parameters are the ones required for my tests (and no others) and that my tests run as before.

Cheers.

[jmichel80]

Humm @fjclark could you investigate the test failure ? We don't want users to inadvertently use geometric combining rules if those have been broken

[fjclark]

@jmichel80 of course.

[lohedges]

Note that I've just updated the test script shown in my comment above. A few of the parameters that I needed to edit in test_combRules.py weren't pasted in the original version for some reason, e.g. turn off thermostat and constraint.

[lohedges]

Okay, I've figured it out, and the current version of SOMD is indeed broken. As @jmichel80 alluded to, combining rules are indeed not being set in SOMD. Comparing with the hard-coded version of OpenMMMD.py in the test script I find the following omission from devel. (Potentially others, who knows.)

git diff                                                                                                            devel ✭ ✱ ◼
diff --git a/wrapper/Tools/OpenMMMD.py b/wrapper/Tools/OpenMMMD.py
index dfd645eb..2117a93f 100644
--- a/wrapper/Tools/OpenMMMD.py
+++ b/wrapper/Tools/OpenMMMD.py
@@ -1208,6 +1208,7 @@ def setupMovesFreeEnergy(system, debug_seed, GPUS, lam_val):
     Integrator_OpenMM.setIntegrator(integrator_type.val)
     Integrator_OpenMM.setFriction(inverse_friction.val)  # Only meaningful for Langevin/Brownian integrators
     Integrator_OpenMM.setPlatform(platform.val)
+    Integrator_OpenMM.setCombiningRules(combining_rules.val)
     Integrator_OpenMM.setConstraintType(constraint.val)
     Integrator_OpenMM.setCutoffType(cutoff_type.val)
     Integrator_OpenMM.setFieldDielectric(rf_dielectric.val)

I can add this to devel now. Would the same thing need to be added to the regular setupMoves, i.e. for somd, not somd-freenrg too, or is this not supported.

[jmichel80]

Thanks for tracking the issue. Edits done OpenMMMD.py by Sofia must have been overwritten at some point. There is no support for geometric combining rules in vanilla somd currently.

[lohedges]

If I git blame that line it was never added in the first place, so I don't think it was overwritten. Looking at the PR here it seems that only the C++ code was updated, not OpenMMD.py.

I'll fix in devel now.

[lohedges]

Yes, it looks like the edits were made in an earlier PR, then reverted following some breaking changes to SOMD. The correct C++ changes were reimplemented in the PR I link to above, but OpenMMMD.py was never restored. It looks like that was the only thing that needed to be re-applied, so things should hopefully be fixed.

[fjclark]

@lohedges can you please confirm that the script you posted is completely updated (as my version fails for the arithmetic test)? Thanks.

[lohedges]

Hmmm, same here. Will check what parameter is missing and update. Apologies.

[lohedges]

Okay, easy fix, not all parameters were being updated in the OpenMMMD namespace. I've updated now. (Must have fixed locally and forgotten to update the code snippet.)

[fjclark]

Hi, I would it be possible to merge this in now that the Boresch restraints PR has been merged into the Sire devel? I've updated it to be consistent with recent changes to Sire and I've checked that it passes. Thanks.

[lohedges]

Hi @fjclark. I don't think this test suite is currently being run. In fact, I don't think the devel branch here works, so this should ideall be merged with the sire-legacy that I created. I'll let @chryswoods comment as to how tests should be managed going forward? Perhaps this would be better placed in this directory?

[chryswoods]

Thanks - yes, I'll merge your changes for sire into the openbiosim/sire repo, and will then test that against the tests here.

I will look at the test and see how to move it into the new pytests test directory. I am going to keep these SireUnitTests alive, as they are useful to run against sire just to make sure we are not badly breaking old code. For now, feel free to keep using this. I'll keep you informed as I move things to openbiosim/sire :-)

[fjclark]

Great, thanks very much both.