Added temperature adjustment to MC barostat for simulated tempering (#3614)

* Added temperature adjustment to MC barostat for simulated tempering

fixes #3612

* Added test for simulated tempering with MC barostat.

* Fixed typo in TestSimulatedTempering.py

* Updated NPT test for SimulatedTempering

* Updated variable naming

* Reverted back to harmonic oscillator test

wrappers/python/openmm/app/simulatedtempering.py

@@ -160,6 +160,9 @@ class SimulatedTempering(object):
         
         self.currentTemperature = 0
         self.simulation.integrator.setTemperature(self.temperatures[self.currentTemperature])
+        for param in self.simulation.context.getParameters():
+            if 'MonteCarloTemperature' in param:
+                self.simulation.context.setParameter(param, self.temperatures[self.currentTemperature])
         
         # Add a reporter to the simulation which will handle the updates and reports.
         
@@ -230,6 +233,9 @@ class SimulatedTempering(object):
                     self._hasMadeTransition = True
                     self.currentTemperature = j
                     self.simulation.integrator.setTemperature(self.temperatures[j])
+                    for param in self.simulation.context.getParameters():
+                        if 'MonteCarloTemperature' in param:
+                            self.simulation.context.setParameter(param, self.temperatures[self.currentTemperature])
                 if self._updateWeights:
                     # Update the weight factors.
                     

wrappers/python/tests/TestSimulatedTempering.py

@@ -53,3 +53,48 @@ class TestSimulatedTempering(unittest.TestCase):
             meanEnergy = sum([0.5*1000*(r-1)**2 for r in d])/n
             expectedEnergy = (0.5*MOLAR_GAS_CONSTANT_R*t).value_in_unit(kilojoules_per_mole)
             self.assertAlmostEqual(expectedEnergy, meanEnergy, delta=expectedEnergy*0.3)
+
+
+    def testHarmonicOscillatorNPT(self):
+        """Test running simulated tempering on a harmonic oscillator with a Monte Carlo barostat."""
+        system = System()
+        system.addParticle(1.0)
+        system.addParticle(1.0)
+        force = HarmonicBondForce()
+        force.addBond(0, 1, 1.0, 1000.0)
+        system.addForce(force)
+        mcbarostat = MonteCarloBarostat(1*bar, 100*kelvin, 2)
+        system.addForce(mcbarostat)
+        force.setUsesPeriodicBoundaryConditions(True)
+        integrator = LangevinIntegrator(100*kelvin, 10/picosecond, 0.001*picosecond)
+        topology = Topology()
+        chain = topology.addChain()
+        residue = topology.addResidue('H2', chain)
+        topology.addAtom('H1', element.hydrogen, residue)
+        topology.addAtom('H2', element.hydrogen, residue)
+        simulation = Simulation(topology, system, integrator, Platform.getPlatformByName('Reference'))
+        simulation.context.setPositions([Vec3(0, 0, 0), Vec3(1, 0, 0)])
+
+        # Check the temperature is correct before creating the simulated tempering simulation
+        
+        self.assertEqual(100*kelvin, integrator.getTemperature())
+        self.assertEqual(100*kelvin, simulation.context.getParameter('MonteCarloTemperature')*kelvin)
+
+        st = SimulatedTempering(simulation, numTemperatures=10, minTemperature=200*kelvin, maxTemperature=400*kelvin, tempChangeInterval=4, reportInterval=10000)
+
+        # Check the temperatures of the integrator and barostat are set to the value of minTemperature
+
+        self.assertEqual(10, len(st.temperatures))
+        self.assertEqual(200*kelvin, st.temperatures[0])
+        self.assertEqual(400*kelvin, st.temperatures[-1])
+        self.assertEqual(200*kelvin, integrator.getTemperature())
+        self.assertEqual(200*kelvin, simulation.context.getParameter('MonteCarloTemperature')*kelvin)
+
+        # Let the simulation run and assert at every step T(mcbarostat) == T(integrator) == T(tempering)
+
+        for i in range(1000):
+            st.step(2)
+            self.assertEqual(st.temperatures[st.currentTemperature], integrator.getTemperature())
+            self.assertEqual(st.temperatures[st.currentTemperature], simulation.context.getParameter('MonteCarloTemperature')*kelvin)
+
+