Merge pull request #264 from choderalab/ideal-gas

Fix ideal gas failures

devtools/conda-recipe/meta.yaml

@@ -19,7 +19,7 @@ requirements:
     - openmm
     - parmed
     - mdtraj
-
+    
   run:
     - python
     - cython

openmmtools/mcmc.py

@@ -288,7 +288,7 @@ def minimize(self, tolerance=1.0*unit.kilocalories_per_mole/unit.angstroms,
         # Retrieve data.
         self.sampler_state.update_from_context(context)
 
-        timer.report_timing()
+        #timer.report_timing()
 
 
 # =============================================================================
@@ -656,8 +656,8 @@ def apply(self, thermodynamic_state, sampler_state):
         timer.start("{}: Context request".format(move_name))
         context, integrator = context_cache.get_context(thermodynamic_state, integrator)
         timer.stop("{}: Context request".format(move_name))
-        logger.debug("{}: Context obtained, platform is {}".format(
-            move_name, context.getPlatform().getName()))
+        #logger.debug("{}: Context obtained, platform is {}".format(
+        #    move_name, context.getPlatform().getName()))
 
         # Perform the integration.
         for attempt_counter in range(self.n_restart_attempts + 1):
@@ -714,7 +714,7 @@ def apply(self, thermodynamic_state, sampler_state):
         sampler_state.update_from_context(context_state)
         timer.stop("{}: update sampler state".format(move_name))
 
-        timer.report_timing()
+        #timer.report_timing()
 
     @abc.abstractmethod
     def _get_integrator(self, thermodynamic_state):
@@ -891,7 +891,7 @@ def apply(self, thermodynamic_state, sampler_state):
 
         # Print timing information.
         timer.stop(benchmark_id)
-        timer.report_timing()
+        #timer.report_timing()
 
     def __getstate__(self):
         if self.context_cache is None:

openmmtools/tests/test_mcmc.py

@@ -91,8 +91,7 @@ def subtest_mcmc_expectation(testsystem, move):
 
     # Test settings.
     temperature = 298.0 * unit.kelvin
-    nequil = 10  # number of equilibration iterations
-    niterations = 40  # number of production iterations
+    niterations = 500  # number of production iterations
     if system.usesPeriodicBoundaryConditions():
         pressure = 1.0*unit.atmosphere
     else:
@@ -109,9 +108,8 @@ def subtest_mcmc_expectation(testsystem, move):
                                              temperature=temperature,
                                              pressure=pressure)
 
-    # Create MCMC sampler and equilibrate.
+    # Create MCMC sampler
     sampler = MCMCSampler(thermodynamic_state, sampler_state, move=move)
-    sampler.run(nequil)
 
     # Accumulate statistics.
     x_n = np.zeros([niterations], np.float64)  # x_n[i] is the x position of atom 1 after iteration i, in angstroms
@@ -143,17 +141,17 @@ def subtest_mcmc_expectation(testsystem, move):
             testsystem.__class__.__name__)
 
         potential_expectation = testsystem.get_potential_expectation(thermodynamic_state) / kT
-        potential_mean = potential_n.mean()
-        g = timeseries.statisticalInefficiency(potential_n, fast=True)
-        dpotential_mean = potential_n.std() / np.sqrt(niterations / g)
+        [t0, g, Neff_max] = timeseries.detectEquilibration(potential_n)
+        potential_mean = potential_n[t0:].mean()
+        dpotential_mean = potential_n[t0:].std() / np.sqrt(Neff_max)
         potential_error = potential_mean - potential_expectation
         nsigma = abs(potential_error) / dpotential_mean
 
         err_msg = ('Potential energy expectation\n'
                    'observed {:10.5f} +- {:10.5f}kT | expected {:10.5f} | '
-                   'error {:10.5f} +- {:10.5f} ({:.1f} sigma)\n'
+                   'error {:10.5f} +- {:10.5f} ({:.1f} sigma) | t0 {:5d} | g {:5.1f} | Neff {:8.1f}\n'
                    '----------------------------------------------------------------------------').format(
-            potential_mean, dpotential_mean, potential_expectation, potential_error, dpotential_mean, nsigma)
+            potential_mean, dpotential_mean, potential_expectation, potential_error, dpotential_mean, nsigma, t0, g, Neff_max)
         assert nsigma <= NSIGMA_CUTOFF, err_msg.format()
         if debug:
             print(err_msg)
@@ -166,17 +164,17 @@ def subtest_mcmc_expectation(testsystem, move):
             testsystem.__class__.__name__)
 
         volume_expectation = testsystem.get_volume_expectation(thermodynamic_state) / (unit.nanometers**3)
-        volume_mean = volume_n.mean()
-        g = timeseries.statisticalInefficiency(volume_n, fast=True)
-        dvolume_mean = volume_n.std() / np.sqrt(niterations / g)
+        [t0, g, Neff_max] = timeseries.detectEquilibration(volume_n)
+        volume_mean = volume_n[t0:].mean()
+        dvolume_mean = volume_n[t0:].std() / np.sqrt(Neff_max)
         volume_error = volume_mean - volume_expectation
         nsigma = abs(volume_error) / dvolume_mean
 
         err_msg = ('Volume expectation\n'
                    'observed {:10.5f} +- {:10.5f}kT | expected {:10.5f} | '
-                   'error {:10.5f} +- {:10.5f} ({:.1f} sigma)\n'
+                   'error {:10.5f} +- {:10.5f} ({:.1f} sigma) | t0 {:5d} | g {:5.1f} | Neff {:8.1f}\n'
                    '----------------------------------------------------------------------------').format(
-            volume_mean, dvolume_mean, volume_expectation, volume_error, dvolume_mean, nsigma)
+            volume_mean, dvolume_mean, volume_expectation, volume_error, dvolume_mean, nsigma, t0, g, Neff_max)
         assert nsigma <= NSIGMA_CUTOFF, err_msg.format()
         if debug:
             print(err_msg)

openmmtools/testsystems.py

@@ -658,7 +658,7 @@ class HarmonicOscillator(TestSystem):
     Notes
     -----
 
-    The natural period of a harmonic oscillator is T = sqrt(m/K), so you will want to use an
+    The natural period of a harmonic oscillator is T = 2*pi*sqrt(m/K), so you will want to use an
     integration timestep smaller than ~ T/10.
 
     The standard deviation in position in each dimension is sigma = (kT / K)^(1/2)