Add support for OpenMM MonteCarloMembraneBarostat

doc/source/changelog.rst

@@ -29,6 +29,8 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * 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 <https://github.com/openbiosim/sire/compare/2025.1.0...2025.2.0>`__ - October 2025
 --------------------------------------------------------------------------------------------
 

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.                 |

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)
 

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

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

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)