From 7a2cb418a7e237dd2c92aea0d24cc33ec51f6984 Mon Sep 17 00:00:00 2001 From: Christopher Woods Date: Fri, 1 Dec 2023 19:59:47 +0000 Subject: [PATCH 1/2] I've made sure that sigma is not set to 0 for any LJ parameters or exceptions. I've added some unit tests for this based on serialising the OpenMM context to XML. I have also checked that the `include_constrained_energies` option works, and added documentation about this plus a unit test. --- doc/source/changelog.rst | 18 +- doc/source/cheatsheet/openmm.rst | 177 +++++++++--------- tests/conftest.py | 21 +-- tests/convert/test_openmm.py | 63 ++++++- wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 19 +- 5 files changed, 187 insertions(+), 111 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 668d2f33e..56a20a7ca 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -19,15 +19,15 @@ organisation on `GitHub `__. a water topology swap. * Added a new :mod:`sire.options` module that contains new - :cls:`sire.options.Option` objects to represent configurable options. + :class:`sire.options.Option` objects to represent configurable options. These include documentation, and make it easier to validate and expose possible values of configurable options. The API docs for - :cls:`~sire.options.Option` shows how to create your own Option type. + :class:`~sire.options.Option` shows how to create your own Option type. The unit test in ``tests/options/test_options.py`` show how to use the options. This is integrated into the sire/OpenMM layer. -* Have :cls:`~sire.io.parser.RST7` return a list of angles from the - ``box_angles()`` function, rather than a :cls:`~sire.maths.Vector`. +* Have :class:`~sire.io.parser.RST7` return a list of angles from the + ``box_angles()`` function, rather than a :class:`~sire.maths.Vector`. This prevents the confusing behaviour where the angles are wrongly shown in units of angstroms... This fixes issues #106. @@ -51,7 +51,11 @@ organisation on `GitHub `__. * Fix validation of ``perturbable_constraint`` dynamics option when the string includes hyphens. This fixes issue #130. -* Fix streaming of :cls:`~sire.vol.TriclinicBox` objects. This fixes issue #128. +* Fix streaming of :class:`~sire.vol.TriclinicBox` objects. This fixes issue #128. + +* Fix the sire to OpenMM conversion so that null LJ parameters will never have + a zero sigma value. They will either be sigma=1/epsilon=0 for non-perturbable + atoms, or sigma=1e-9/epsilon=1e-9 for perturbable atoms. * Please add an item to this changelog when you create your PR @@ -145,7 +149,7 @@ organisation on `GitHub `__. * Added support for alchemical restraints to the OpenMM dynamics layer. This lets you scale restraints as part of a λ-coordinate. This is - documented in the :doc:`tutorial `. Restraints can be named, meaning that you can scale different restraints at different stages and by different values across the λ-coordinate. @@ -448,7 +452,7 @@ organisation on `GitHub `__. * Added support for performing minimisation and molecular dynamics simulations based on the OpenMM integration. This is documented in full via both - :doc:`a tutorial ` and a + :doc:`a tutorial ` and a :doc:`detailed guide `. * Fixed the Amber PRMTOP `dihedral ring bug `__. diff --git a/doc/source/cheatsheet/openmm.rst b/doc/source/cheatsheet/openmm.rst index d5345681c..40cb191ec 100644 --- a/doc/source/cheatsheet/openmm.rst +++ b/doc/source/cheatsheet/openmm.rst @@ -73,91 +73,97 @@ by setting values for the ``temperature`` and ``pressure`` keys. Available keys and allowable values are listed below. -+---------------------------+----------------------------------------------------------+ -| Key | Valid values | -+===========================+==========================================================+ -| com_reset_frequency | The frequency at which the ``CMMotionRemover`` acts to | -| | remove center of mass relative motion. If this is not | -| | set (the default) then center of mass motion is not | -| | removed. | -+---------------------------+----------------------------------------------------------+ -| constraint | Type of constraint to use for bonds and/or angles. | -| | Valid strings are ``none``, ``h-bonds``, ``bonds``, | -| | ``h-bonds-h-angles`` and ``bonds-h-angles``. | -+---------------------------+----------------------------------------------------------+ -| coulomb_power | The coulomb power parameter used by the softening | -| | potential used to soften interactions involving | -| | ghost atoms. | -+---------------------------+----------------------------------------------------------+ -| cutoff | Size of the non-bonded cutoff, e.g. | -| | ``7.5*sr.units.angstrom`` | -+---------------------------+----------------------------------------------------------+ -| cutoff_type | Type of cutoff, e.g. ``PARTICLE_MESH_EWALD``, ``PME``, | -| | ``NO_CUTOFF``, ``REACTION_FIELD``, ``RF``, ``EWALD`` | -+---------------------------+----------------------------------------------------------+ -| cpu_pme | Boolean value, e.g. ``True`` or ``False`` as to whether | -| | or not PME should be evaluated on the CPU rather than | -| | on the GPU. | -+---------------------------+----------------------------------------------------------+ -| device | Any valid OpenMM device number or device string, e.g. | -| | ``0`` would select device 0 | -+---------------------------+----------------------------------------------------------+ -| dielectric | Dielectric value if a reaction field cutoff is used, | -| | e.g. ``78.0`` | -+---------------------------+----------------------------------------------------------+ -| fixed | The atoms in the system that should be fixed (not moved) | -+---------------------------+----------------------------------------------------------+ -| ignore_perturbations | Whether or not to ignore any perturbations and only set | -| | up a perturbable molecule as a non-perurbable molecule | -| | from only the reference state. | -+---------------------------+----------------------------------------------------------+ -| integrator | The MD integrator to use, e.g. | -| | ``verlet``, ``leapfrog``, ``langevin``, | -| | ``langevin_middle``, ``nose_hoover``, | -| | ``brownian``, ``andersen`` | -+---------------------------+----------------------------------------------------------+ -| friction | Friction value for the integrator, in inverse time, e.g. | -| | ``5.0 / sr.units.picosecond`` | -+---------------------------+----------------------------------------------------------+ -| lambda | The λ-value at which to set up the system (assuming this | -| | contains any perturbable molecules or restraints) | -+---------------------------+----------------------------------------------------------+ -| platform | Any valid OpenMM platform string, e.g. ``CUDA``, | -| | ``OpenCL``, ``Metal``, ```CPU``, ``Reference`` | -+---------------------------+----------------------------------------------------------+ -| precision | Any valid OpenMM platform precision value, e.g. | -| | ``single``, ``mixed`` or ``double``. | -+---------------------------+----------------------------------------------------------+ -| pressure | Any pressure value, e.g. ``1*sr.units.atm`` | -+---------------------------+----------------------------------------------------------+ -| restraints | The :class:`~sire.mm.Restraints` object (or list of | -| | objects) that describe the restraints that should be | -| | added to the system. | -+---------------------------+----------------------------------------------------------+ -| schedule | The :class:`~sire.cas.LambdaSchedule` to use that | -| | controls how parameters are modified with λ | -+---------------------------+----------------------------------------------------------+ -| shift_delta | The shift_delta parameter to use for the softening | -| | potential used to soften interactions involving | -| | ghost atoms. | -+---------------------------+----------------------------------------------------------+ -| space | Space in which the simulation should be conducted, e.g. | -| | `sr.vol.Cartesian` | -+---------------------------+----------------------------------------------------------+ -| swap_end_states | Whether to swap the end states of a perturbable molecule | -| | (i.e. treat the perturbed state as the reference state | -| | and vice versa). | -+---------------------------+----------------------------------------------------------+ -| temperature | Any temperature value, e.g. ``25*sr.units.celsius`` | -+---------------------------+----------------------------------------------------------+ -| threads | The number of threads to use in the CPU platform | -+---------------------------+----------------------------------------------------------+ -| timestep | Time between integration steps, e.g. | -| | ``2 * sr.units.femtosecond`` | -+---------------------------+----------------------------------------------------------+ -| use_dispersion_correction | Whether or not to use the dispersion correction to | -| | deal with cutoff issues. This is very expensive. | -+---------------------------+----------------------------------------------------------+ ++------------------------------+----------------------------------------------------------+ +| Key | Valid values | ++==============================+==========================================================+ +| com_reset_frequency | The frequency at which the ``CMMotionRemover`` acts to | +| | remove center of mass relative motion. If this is not | +| | set (the default) then center of mass motion is not | +| | removed. | ++------------------------------+----------------------------------------------------------+ +| constraint | Type of constraint to use for bonds and/or angles. | +| | Valid strings are ``none``, ``h-bonds``, ``bonds``, | +| | ``h-bonds-h-angles`` and ``bonds-h-angles``. | ++------------------------------+----------------------------------------------------------+ +| coulomb_power | The coulomb power parameter used by the softening | +| | potential used to soften interactions involving | +| | ghost atoms. | ++------------------------------+----------------------------------------------------------+ +| cutoff | Size of the non-bonded cutoff, e.g. | +| | ``7.5*sr.units.angstrom`` | ++------------------------------+----------------------------------------------------------+ +| cutoff_type | Type of cutoff, e.g. ``PARTICLE_MESH_EWALD``, ``PME``, | +| | ``NO_CUTOFF``, ``REACTION_FIELD``, ``RF``, ``EWALD`` | ++------------------------------+----------------------------------------------------------+ +| cpu_pme | Boolean value, e.g. ``True`` or ``False`` as to whether | +| | or not PME should be evaluated on the CPU rather than | +| | on the GPU. | ++------------------------------+----------------------------------------------------------+ +| device | Any valid OpenMM device number or device string, e.g. | +| | ``0`` would select device 0 | ++------------------------------+----------------------------------------------------------+ +| dielectric | Dielectric value if a reaction field cutoff is used, | +| | e.g. ``78.0`` | ++------------------------------+----------------------------------------------------------+ +| fixed | The atoms in the system that should be fixed (not moved) | ++------------------------------+----------------------------------------------------------+ +| ignore_perturbations | Whether or not to ignore any perturbations and only set | +| | up a perturbable molecule as a non-perurbable molecule | +| | from only the reference state. | ++------------------------------+----------------------------------------------------------+ +| include_constrained_energies | Whether or not to include the bond and angle energies | +| | of bonds and angles that are constrained. | +| | This defaults to ``True``, as we normally do want | +| | to calculate all of the energies of the internals, | +| | even if they are constrained. | ++------------------------------+----------------------------------------------------------+ +| integrator | The MD integrator to use, e.g. | +| | ``verlet``, ``leapfrog``, ``langevin``, | +| | ``langevin_middle``, ``nose_hoover``, | +| | ``brownian``, ``andersen`` | ++------------------------------+----------------------------------------------------------+ +| friction | Friction value for the integrator, in inverse time, e.g. | +| | ``5.0 / sr.units.picosecond`` | ++------------------------------+----------------------------------------------------------+ +| lambda | The λ-value at which to set up the system (assuming this | +| | contains any perturbable molecules or restraints) | ++------------------------------+----------------------------------------------------------+ +| platform | Any valid OpenMM platform string, e.g. ``CUDA``, | +| | ``OpenCL``, ``Metal``, ```CPU``, ``Reference`` | ++------------------------------+----------------------------------------------------------+ +| precision | Any valid OpenMM platform precision value, e.g. | +| | ``single``, ``mixed`` or ``double``. | ++------------------------------+----------------------------------------------------------+ +| pressure | Any pressure value, e.g. ``1*sr.units.atm`` | ++------------------------------+----------------------------------------------------------+ +| restraints | The :class:`~sire.mm.Restraints` object (or list of | +| | objects) that describe the restraints that should be | +| | added to the system. | ++------------------------------+----------------------------------------------------------+ +| schedule | The :class:`~sire.cas.LambdaSchedule` to use that | +| | controls how parameters are modified with λ | ++------------------------------+----------------------------------------------------------+ +| shift_delta | The shift_delta parameter to use for the softening | +| | potential used to soften interactions involving | +| | ghost atoms. | ++------------------------------+----------------------------------------------------------+ +| space | Space in which the simulation should be conducted, e.g. | +| | `sr.vol.Cartesian` | ++------------------------------+----------------------------------------------------------+ +| swap_end_states | Whether to swap the end states of a perturbable molecule | +| | (i.e. treat the perturbed state as the reference state | +| | and vice versa). | ++------------------------------+----------------------------------------------------------+ +| temperature | Any temperature value, e.g. ``25*sr.units.celsius`` | ++------------------------------+----------------------------------------------------------+ +| threads | The number of threads to use in the CPU platform | ++------------------------------+----------------------------------------------------------+ +| timestep | Time between integration steps, e.g. | +| | ``2 * sr.units.femtosecond`` | ++------------------------------+----------------------------------------------------------+ +| use_dispersion_correction | Whether or not to use the dispersion correction to | +| | deal with cutoff issues. This is very expensive. | ++------------------------------+----------------------------------------------------------+ Higher level API ---------------- @@ -336,4 +342,3 @@ from files (e.g. via that property of the ``.trajectory()`` function). The energies are saved to the ``energy_trajectory`` property of the returned molecules, and accessible via that property or the :func:`~sire.system.System.energy_trajectory` function. - diff --git a/tests/conftest.py b/tests/conftest.py index 88f267b8a..99409d09b 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -43,9 +43,7 @@ def pytest_collection_modifyitems(config, items): item.add_marker(skip_slow) if not runveryslow: - skip_veryslow = pytest.mark.skip( - reason="need --runveryslow option to run" - ) + skip_veryslow = pytest.mark.skip(reason="need --runveryslow option to run") for item in items: if "veryslow" in item.keywords: item.add_marker(skip_veryslow) @@ -113,23 +111,17 @@ def wrapped_mols(): @pytest.fixture(scope="session") def openmm_interchange_mols(): - return sr.load_test_files( - "openmm_interchange.rst7", "openmm_interchange.prm7" - ) + return sr.load_test_files("openmm_interchange.rst7", "openmm_interchange.prm7") @pytest.fixture(scope="session") def triclinic_protein(): - return sr.load_test_files( - "triclinic_protein.prm7", "triclinic_protein.crd" - ) + return sr.load_test_files("triclinic_protein.prm7", "triclinic_protein.crd") @pytest.fixture(scope="session") def triclinic_protein_rst7(): - return sr.load_test_files( - "triclinic_protein.prm7", "triclinic_protein.rst" - ) + return sr.load_test_files("triclinic_protein.prm7", "triclinic_protein.rst") @pytest.fixture(scope="session") @@ -150,3 +142,8 @@ def ethane_12dichloroethane(): @pytest.fixture(scope="session") def pentane_cyclopentane(): return sr.load_test_files("pentane_cyclopentane.bss") + + +@pytest.fixture(scope="session") +def zero_lj_mols(): + return sr.load_test_files("zero_lj.prm7", "zero_lj.rst7") diff --git a/tests/convert/test_openmm.py b/tests/convert/test_openmm.py index de104ecce..a34757c72 100644 --- a/tests/convert/test_openmm.py +++ b/tests/convert/test_openmm.py @@ -39,9 +39,7 @@ def test_openmm_single_energy(kigaki_mols): energy = energy.value_in_unit(energy.unit) # these won't be exactly the same - this is 5227 +/- 0.1 kJ mol-1 - assert mol.energy(map=map).to(sr.units.kJ_per_mol) == pytest.approx( - energy, abs=0.1 - ) + assert mol.energy(map=map).to(sr.units.kJ_per_mol) == pytest.approx(energy, abs=0.1) @pytest.mark.skipif( @@ -285,3 +283,62 @@ def test_openmm_ignore_constrained(ala_mols): # these two energies should be different, because # we should be ignoring the constrained bonds and angles assert abs(nrg2.value() - nrg1.value()) > 1.0 + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_no_zero_sigmas(zero_lj_mols): + mols = zero_lj_mols + + omm = sr.convert.to(mols, "openmm", map={"constraint": "h-bonds"}) + + from openmm import XmlSerializer + + xml = XmlSerializer.serialize(omm.getSystem()) + + for line in xml.split("\n"): + assert 'sig="0"' not in line + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_skipped_constrained_bonds(zero_lj_mols): + mols = zero_lj_mols + + omm1 = sr.convert.to( + mols, + "openmm", + map={"constraint": "h-bonds", "include_constrained_energies": True}, + ) + + omm2 = sr.convert.to( + mols, + "openmm", + map={"constraint": "h-bonds", "include_constrained_energies": False}, + ) + + from openmm import XmlSerializer + + xml1 = XmlSerializer.serialize(omm1.getSystem()) + xml2 = XmlSerializer.serialize(omm2.getSystem()) + + assert xml1 != xml2 + + lines1 = xml1.split("\n") + lines2 = xml2.split("\n") + + i = 0 + + for j in range(0, len(lines2)): + line1 = lines1[i] + line2 = lines2[j] + i += 1 + + while line1 != line2: + line1 = lines1[i] + assert "Bond" in line1 + i += 1 diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 048a62c5d..dbc1f3276 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -296,7 +296,7 @@ std::tuple OpenMMMolecule::getException( epsilon = lj_14_scl * std::sqrt(std::get<2>(clj0) * std::get<2>(clj1)); } - if (this->isPerturbable() and charge == 0 and epsilon == 0) + if (this->isPerturbable() and charge == 0 and std::abs(epsilon) < 1e-9) { // openmm tries to optimise away zero parameters - this is an issue // as perturbation requires that we don't remove them! @@ -306,6 +306,12 @@ std::tuple OpenMMMolecule::getException( sigma = 1e-9; epsilon = 1e-9; } + else if (sigma == 0) + { + // make sure we never have zero sigma values - this is a null parameter + sigma = 1; + epsilon = 0; + } return std::make_tuple(atom0 + start_index, atom1 + start_index, @@ -491,8 +497,15 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, const double chg = params_charges.at(idx_to_cgatomidx_data[i]).to(SireUnits::mod_electron); const auto &lj = params_ljs.at(idx_to_cgatomidx_data[i]); - const double sig = lj.sigma().to(SireUnits::nanometer); - const double eps = lj.epsilon().to(SireUnits::kJ_per_mol); + double sig = lj.sigma().to(SireUnits::nanometer); + double eps = lj.epsilon().to(SireUnits::kJ_per_mol); + + if (std::abs(sig) < 1e-9) + { + // this must be a null parameter + eps = 0; + sig = 1; + } cljs_data[i] = std::make_tuple(chg, sig, eps); } From 52d3b8765791399e9df093c3df6bb0ce9c6e409a Mon Sep 17 00:00:00 2001 From: Christopher Woods Date: Sat, 2 Dec 2023 10:55:25 +0000 Subject: [PATCH 2/2] Fixed platform issue that caused test failure on GH MacOS runner Added in an energy check so that we test for regressions in calculating the total energy for contexts with and without energy terms for constraint internals --- tests/convert/test_openmm.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/tests/convert/test_openmm.py b/tests/convert/test_openmm.py index a34757c72..4a014ebf5 100644 --- a/tests/convert/test_openmm.py +++ b/tests/convert/test_openmm.py @@ -292,7 +292,9 @@ def test_openmm_ignore_constrained(ala_mols): def test_openmm_no_zero_sigmas(zero_lj_mols): mols = zero_lj_mols - omm = sr.convert.to(mols, "openmm", map={"constraint": "h-bonds"}) + omm = sr.convert.to(mols, "openmm", + map={"constraint": "h-bonds", + "platform": "Reference"}) from openmm import XmlSerializer @@ -312,15 +314,27 @@ def test_openmm_skipped_constrained_bonds(zero_lj_mols): omm1 = sr.convert.to( mols, "openmm", - map={"constraint": "h-bonds", "include_constrained_energies": True}, + map={"constraint": "h-bonds", + "include_constrained_energies": True, + "platform": "Reference"}, ) omm2 = sr.convert.to( mols, "openmm", - map={"constraint": "h-bonds", "include_constrained_energies": False}, + map={"constraint": "h-bonds", + "include_constrained_energies": False, + "platform": "Reference"}, ) + nrg1 = omm1.get_potential_energy().to(sr.units.kcal_per_mol) + nrg2 = omm2.get_potential_energy().to(sr.units.kcal_per_mol) + + # Check the energies haven't changed + # (regression check - here are the current values) + assert nrg1 == pytest.approx(-447.44, 1e-3) + assert nrg2 == pytest.approx(-3279.87, 1e-3) + from openmm import XmlSerializer xml1 = XmlSerializer.serialize(omm1.getSystem())