From 032b4ddecca8fc167efa29ab7db81df2bc05f138 Mon Sep 17 00:00:00 2001 From: Christopher Woods Date: Tue, 5 Dec 2023 17:23:52 +0000 Subject: [PATCH] Backported the fix in PR 132 to main --- doc/source/cheatsheet/openmm.rst | 177 +++++++++--------- tests/conftest.py | 5 + tests/convert/test_openmm.py | 73 ++++++++ wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 19 +- 4 files changed, 185 insertions(+), 89 deletions(-) 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..06adf0b75 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -150,3 +150,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..ed04ac302 100644 --- a/tests/convert/test_openmm.py +++ b/tests/convert/test_openmm.py @@ -285,3 +285,76 @@ 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", + "platform": "Reference"}) + + 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, + "platform": "Reference"}, + ) + + omm2 = sr.convert.to( + mols, + "openmm", + 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()) + 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 9df2420a9..d21c03ae7 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); }