Merge pull request #1049 from markovmodel/thermo-doc

Thermo doc

devtools/conda-recipe/meta.yaml

@@ -33,7 +33,7 @@ requirements:
     - scipy
     - setuptools
     - six >=1.10
-    - thermotools >=0.2.3
+    - thermotools >=0.2.5
 
 test:
   requires:

doc/source/CHANGELOG.rst

@@ -3,6 +3,13 @@ Changelog
 
 2.3.2 (tba)
 
+**New features**:
+
+- thermo:
+   - Allow for periodicity in estimate_umbrella_sampling()
+   - Add *_full_state getter variants to access stationary properties on the full set of states
+     instead of the active set
+
 **Fixes**:
 
 - coordinates:

pyemma/thermo/api.py

@@ -43,7 +43,7 @@ def estimate_umbrella_sampling(
     us_trajs, us_dtrajs, us_centers, us_force_constants, md_trajs=None, md_dtrajs=None, kT=None,
     maxiter=10000, maxerr=1.0E-15, save_convergence_info=0,
     estimator='wham', lag=1, dt_traj='1 step', init=None, init_maxiter=10000, init_maxerr=1.0E-8,
-    **kwargs):
+    width=None, **kwargs):
     r"""
     This function acts as a wrapper for ``tram()``, ``dtram()``, ``mbar``, and ``wham()`` and
     handles the calculation of bias energies (``bias``) and thermodynamic state trajectories
@@ -112,6 +112,10 @@ def estimate_umbrella_sampling(
         The maximum number of self-consistent iterations during the initialization.
     init_maxerr : float, optional, default=1.0E-8
         Convergence criterion for the initialization.
+    width : array-like of float, optional, default=None
+        Specify periodicity for individual us_traj dimensions. Each positive entry will make the
+        corresponding feature periodic and use the given value as width. None/zero values will be
+        treated as non-periodic.
     **kwargs : dict, optional
         You can use this to pass estimator-specific named parameters to the chosen estimator, which
         are not already coverd by ``estimate_umbrella_sampling()``.
@@ -229,7 +233,7 @@ def estimate_umbrella_sampling(
             i += 1
     # data preparation
     ttrajs, btrajs, umbrella_centers, force_constants, unbiased_state = _get_umbrella_sampling_data(
-        us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT)
+        us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT, width=width)
     estimator_obj = None
     # estimation
     if estimator == 'wham':

pyemma/thermo/tests/test_util.py

@@ -317,8 +317,9 @@ def test_umbrella_sampling_bias_sequences_1x0(self):
         force_constants = np.array([
             util._ensure_force_constant(1.0, 1),
             util._ensure_force_constant(2.0, 1)], dtype=np.float64)
+        width = np.zeros(shape=(umbrella_centers.shape[1],), dtype=np.float64)
         self._assert_bias_sequences(
-            util._get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants),
+            util._get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants, width),
             [np.array([[0.0, 1.0], [0.125, 0.25], [0.5, 0.0]])])
 
     def test_umbrella_sampling_bias_sequences_1x1(self):
@@ -329,48 +330,72 @@ def test_umbrella_sampling_bias_sequences_1x1(self):
         force_constants = np.array([
             util._ensure_force_constant(1.0, 1),
             util._ensure_force_constant(2.0, 1)], dtype=np.float64)
+        width = np.zeros(shape=(umbrella_centers.shape[1],), dtype=np.float64)
         self._assert_bias_sequences(
-            util._get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants),
+            util._get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants, width),
             [np.array([[0.0, 1.0], [0.125, 0.25], [0.5, 0.0]])])
 
     def test_umbrella_sampling_bias_sequences_catches_unmatching_dimension(self):
         # wrong centers + constants
         with self.assertRaises(TypeError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                np.array([[1.0, 1.0]]), [[[1.0, 0.0], [1.0, 0.0]]])
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                np.array([[1.0, 1.0]]), [[[1.0, 0.0], [1.0, 0.0]]],
+                np.array([0.0, 0.0]))
         with self.assertRaises(TypeError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                [[1.0, 1.0]], np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                [[1.0, 1.0]], np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))
         with self.assertRaises(ValueError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                np.array([1.0, 1.0]), np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                np.array([1.0, 1.0]), np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))
         with self.assertRaises(ValueError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                np.array([[1.0, 1.0]]), np.array([[1.0, 0.0], [1.0, 0.0]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                np.array([[1.0, 1.0]]), np.array([[1.0, 0.0], [1.0, 0.0]]),
+                np.array([0.0, 0.0]))
         with self.assertRaises(ValueError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                np.array([[[1.0, 1.0]]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                np.array([[[1.0, 1.0]]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))
         with self.assertRaises(ValueError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                np.array([[1.0, 1.0]]), np.array([[[[1.0, 0.0], [1.0, 0.0]]]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                np.array([[1.0, 1.0]]), np.array([[[[1.0, 0.0], [1.0, 0.0]]]]),
+                np.array([0.0, 0.0]))
         # conflicting centers + constants
         with self.assertRaises(ValueError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                np.array([[1.0, 1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                np.array([[1.0, 1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))
         with self.assertRaises(ValueError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                np.array([[1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                np.array([[1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))
         with self.assertRaises(ValueError):
-            util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
-                np.array([[1.0, 1.0], [2.0, 2.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])],
+                np.array([[1.0, 1.0], [2.0, 2.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))
         # traj does not match valid centers + constants
         with self.assertRaises(TypeError):
-            util._get_umbrella_bias_sequences([[[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]]],
-                np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+            util._get_umbrella_bias_sequences(
+                [[[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]]],
+                np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))
         with self.assertRaises(ValueError):
-            util._get_umbrella_bias_sequences([np.array([0.0, 0.5, 1.0])],
-                np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+            util._get_umbrella_bias_sequences(
+                [np.array([0.0, 0.5, 1.0])],
+                np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))
         with self.assertRaises(ValueError):
             util._get_umbrella_bias_sequences(
                 [np.array([[0.0, 1.0, 2.0], [0.5, 1.0, 2.0], [1.0, 1.0, 2.0]])],
-                np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]))
+                np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]]),
+                np.array([0.0, 0.0]))

pyemma/thermo/util/util.py

@@ -185,7 +185,7 @@ def _get_umbrella_sampling_parameters(
         force_constants /= kT
     return ttrajs, umbrella_centers, force_constants, unbiased_state
 
-def _get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants):
+def _get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants, width):
     from thermotools.util import get_umbrella_bias as _get_umbrella_bias
     bias_sequences = []
     if not isinstance(umbrella_centers, _np.ndarray):
@@ -215,10 +215,11 @@ def _get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants):
             raise ValueError("Trajectory %d has unmatching dimension: %d!=%d" % (
                 i, traj.shape[1], dimension))
         bias_sequences.append(
-            _get_umbrella_bias(traj, umbrella_centers, force_constants))
+            _get_umbrella_bias(traj, umbrella_centers, force_constants, width))
     return bias_sequences
 
-def get_umbrella_sampling_data(us_trajs, us_centers, us_force_constants, md_trajs=None, kT=None):
+def get_umbrella_sampling_data(
+    us_trajs, us_centers, us_force_constants, md_trajs=None, kT=None, width=None):
     r"""
     Wraps umbrella sampling data or a mix of umbrella sampling and and direct molecular dynamics.
 
@@ -239,6 +240,10 @@ def get_umbrella_sampling_data(us_trajs, us_centers, us_force_constants, md_traj
         Unbiased molecular dynamics simulations. Format like umbrella_trajs.
     kT : float (optinal)
         Use this attribute if the supplied force constants are NOT unit-less.
+    width : array-like of float, optional, default=None
+        Specify periodicity for individual us_traj dimensions. Each positive entry will make the
+        corresponding feature periodic and use the given value as width. None/zero values will be
+        treated as non-periodic.
 
     Returns
     -------
@@ -258,7 +263,16 @@ def get_umbrella_sampling_data(us_trajs, us_centers, us_force_constants, md_traj
         us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT)
     if md_trajs is None:
         md_trajs = []
-    btrajs = _get_umbrella_bias_sequences(us_trajs + md_trajs, umbrella_centers, force_constants)
+    if width is None:
+        width = _np.zeros(shape=(umbrella_centers.shape[1],), dtype=_np.float64)
+    else:
+        width = _np.asarray(
+            map(lambda w: w if w is not None and w > 0.0 else 0.0, width),
+            dtype=_np.float64)
+    if width.shape[0] != umbrella_centers.shape[1]:
+        raise ValueError('Unmatching number of width components.')
+    btrajs = _get_umbrella_bias_sequences(
+        us_trajs + md_trajs, umbrella_centers, force_constants, width)
     return ttrajs, btrajs, umbrella_centers, force_constants, unbiased_state
 
 # ==================================================================================================

setup.py

@@ -248,7 +248,7 @@ def run_tests(self):
                       'mdtraj>=1.8.0',
                       'matplotlib',
                       'msmtools>=1.2',
-                      'thermotools>=0.2.3',
+                      'thermotools>=0.2.5',
                       'bhmm>=0.6,<0.7',
                       'joblib>0.8.4',
                       'pyyaml',