Fixes missing barostat issue

news/barostat_fix.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* <news item>
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* Fixes missing barostat in non-alchemical simulations for SFE Protocol (issue #114)
+
+**Security:**
+
+* <news item>

src/pontibus/protocols/solvation/base.py

@@ -23,6 +23,7 @@
 from openfe.utils import log_system_probe, without_oechem_backend
 from openff.interchange.interop.openmm import to_openmm_positions
 from openff.toolkit import Molecule as OFFMolecule
+from openff.units.openmm import to_openmm
 from openmm import app
 from openmmtools.alchemy import (
     AbsoluteAlchemicalFactory,
@@ -204,6 +205,15 @@ def _get_omm_objects(
             if isinstance(force, openmm.CMMotionRemover):
                 omm_system.removeForce(idx)
 
+        # Add a barostat if needed
+        if solvent_component is not None:
+            barostat = openmm.MonteCarloBarostat(
+                to_openmm(settings["thermo_settings"].pressure),
+                to_openmm(settings["thermo_settings"].temperature),
+                settings["integrator_settings"].barostat_frequency.m,
+            )
+            omm_system.addForce(barostat)
+
         positions = to_openmm_positions(interchange, include_virtual_sites=True)
 
         # Post creation system validation

src/pontibus/tests/protocols/solvation/test_dry_run.py

@@ -6,7 +6,7 @@
 import pytest
 from gufe import ChemicalSystem
 from openff.units import unit
-from openff.units.openmm import ensure_quantity
+from openff.units.openmm import ensure_quantity, from_openmm
 from openmm import (
     CustomBondForce,
     CustomNonbondedForce,
@@ -90,6 +90,8 @@ def test_dry_run_solv_benzene(experimental, charged_benzene, tmpdir):
     s = ASFEProtocol.default_settings()
     s.protocol_repeats = 1
     s.solvent_output_settings.output_indices = "resname AAA"
+    # Set a random barostat frequency to make sure it goes all the way
+    s.integrator_settings.barostat_frequency = 125
     s.alchemical_settings.experimental = experimental
 
     protocol = ASFEProtocol(
@@ -148,6 +150,13 @@ def assert_force_num(system, forcetype, number):
         assert_force_num(system, PeriodicTorsionForce, 1)
         assert_force_num(system, MonteCarloBarostat, 1)
 
+        # Check the initial barostat made it all the way through
+        for force in system.getForces():
+            if isinstance(force, MonteCarloBarostat):
+                assert force.getFrequency() == 125
+                assert from_openmm(force.getDefaultPressure()) == s.thermo_settings.pressure
+                assert from_openmm(force.getDefaultTemperature()) == s.thermo_settings.temperature
+
         # Check the nonbonded force is PME
         nonbond = [f for f in system.getForces() if isinstance(f, NonbondedForce)]
         assert nonbond[0].getNonbondedMethod() == NonbondedForce.PME