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..4a014ebf5 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,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 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); }