diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 0431accdd..f30a236c3 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -29,6 +29,8 @@ organisation on `GitHub `__. * Added :func:`Sire::IO::setCoordinates` function to set atom coordinates of an entire system. +* Add support for ``openmm.MonteCarloMembraneBarostat`` in Sire-to-OpenMM conversion. + `2025.2.0 `__ - October 2025 -------------------------------------------------------------------------------------------- diff --git a/doc/source/cheatsheet/openmm.rst b/doc/source/cheatsheet/openmm.rst index d53daa436..4f27fccbc 100644 --- a/doc/source/cheatsheet/openmm.rst +++ b/doc/source/cheatsheet/openmm.rst @@ -199,6 +199,9 @@ Available keys and allowable values are listed below. | space | Space in which the simulation should be conducted, e.g. | | | `sr.vol.Cartesian` | +------------------------------+----------------------------------------------------------+ +| surface_tension | Surface tension for membrane simulations at constant | +| | pressure, e.g. ``0.05*sr.units.bar*sr.units.nanometer`` | ++------------------------------+----------------------------------------------------------+ | 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). This defaults to False. | diff --git a/src/sire/mol/__init__.py b/src/sire/mol/__init__.py index f3a7ded7d..2f4d09770 100644 --- a/src/sire/mol/__init__.py +++ b/src/sire/mol/__init__.py @@ -1567,6 +1567,7 @@ def _dynamics( ignore_perturbations=None, temperature=None, pressure=None, + surface_tension=None, rest2_scale=None, rest2_selection=None, vacuum=None, @@ -1688,6 +1689,10 @@ def _dynamics( microcanonical (NVE) or canonical (NVT) simulation will be run if the pressure is not set. + surface_tension: surface_tension + The surface tension to use if running a NPT simulation with + a membrane barostat. + rest2_scale: float The scaling factor to apply when running a REST2 simulation (Replica Exchange with Solute Tempering). This defaults to 1.0. @@ -1878,6 +1883,10 @@ def _dynamics( pressure = u(pressure) map.set("pressure", pressure) + if surface_tension is not None: + surface_tension = u(surface_tension) + map.set("surface_tension", surface_tension) + if rest2_scale is not None: map.set("rest2_scale", rest2_scale) diff --git a/tests/convert/test_openmm.py b/tests/convert/test_openmm.py index 4a6b58d69..d2efb7c15 100644 --- a/tests/convert/test_openmm.py +++ b/tests/convert/test_openmm.py @@ -447,3 +447,45 @@ def test_openmm_default_box_vectors(ala_mols, openmm_platform): # Check that the box vectors match the expected values. for i, vec in enumerate(box_vectors): assert vec[i].value_in_unit(angstrom) == pytest.approx(box[i].value(), abs=1e-3) + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_membrane_barostat(ala_mols, openmm_platform): + from openmm import MonteCarloMembraneBarostat + from openmm import unit + + mols = ala_mols.clone() + + # Create a dynamics object with a membrane barostat. + d = mols.dynamics( + pressure="1atm", + temperature="298K", + surface_tension="1 angstrom*atm", + platform=openmm_platform, + barostat_frequency=50, + ) + + # Find the barostat. + barostat = None + for force in d.context().getSystem().getForces(): + if force.getName() == "MonteCarloMembraneBarostat": + barostat = force + break + assert barostat is not None + + # Check the barostat parameters. + assert barostat.getDefaultPressure().value_in_unit( + unit.atmosphere + ) == pytest.approx(1.0, abs=1e-3) + assert barostat.getDefaultSurfaceTension().value_in_unit( + unit.angstrom * unit.atmosphere + ) == pytest.approx(1.0, abs=1e-3) + assert barostat.getDefaultTemperature().value_in_unit(unit.kelvin) == pytest.approx( + 298.0, abs=1e-3 + ) + assert barostat.getXYMode() == MonteCarloMembraneBarostat.XYIsotropic + assert barostat.getZMode() == MonteCarloMembraneBarostat.ZFree + assert barostat.getFrequency() == 50 diff --git a/wrapper/Convert/SireOpenMM/_sommcontext.py b/wrapper/Convert/SireOpenMM/_sommcontext.py index 15fcba546..8d990cc46 100644 --- a/wrapper/Convert/SireOpenMM/_sommcontext.py +++ b/wrapper/Convert/SireOpenMM/_sommcontext.py @@ -309,6 +309,12 @@ def set_pressure(self, pressure): """ raise NotImplementedError("We can't yet set the pressure") + def set_surface_tension(self, surface_tension): + """ + Set the target surface tension for the dynamics. + """ + raise NotImplementedError("We can't yet set the surface tension") + def get_potential_energy(self, to_sire_units: bool = True): """ Calculate and return the potential energy of the system diff --git a/wrapper/Convert/__init__.py b/wrapper/Convert/__init__.py index a481b4b86..34e706248 100644 --- a/wrapper/Convert/__init__.py +++ b/wrapper/Convert/__init__.py @@ -223,6 +223,12 @@ def sire_to_openmm(mols, map): friction = friction.to(1.0 / picosecond) / openmm.unit.picosecond + if map.specified("surface_tension"): + surface_tension = map["surface_tension"].value() + surface_tension = surface_tension.to("bar * nanometer") + else: + surface_tension = None + use_andersen = False temperature = None @@ -354,7 +360,19 @@ def sire_to_openmm(mols, map): pressure = ensemble.pressure().to(atm) * openmm.unit.atmosphere - barostat = openmm.MonteCarloBarostat(pressure, temperature, barostat_freq) + if surface_tension is not None: + barostat = openmm.MonteCarloMembraneBarostat( + pressure, + surface_tension, + temperature, + openmm.MonteCarloMembraneBarostat.XYIsotropic, + openmm.MonteCarloMembraneBarostat.ZFree, + barostat_freq, + ) + else: + barostat = openmm.MonteCarloBarostat( + pressure, temperature, barostat_freq + ) system.addForce(barostat)