diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml index 7faef6caa..db5738a0e 100644 --- a/devtools/conda-recipe/meta.yaml +++ b/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: diff --git a/doc/source/CHANGELOG.rst b/doc/source/CHANGELOG.rst index 2a1824c68..304c1068c 100644 --- a/doc/source/CHANGELOG.rst +++ b/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: diff --git a/pyemma/thermo/api.py b/pyemma/thermo/api.py index 80c1f3178..bef64f8a7 100644 --- a/pyemma/thermo/api.py +++ b/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': diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index 105e4306c..b67c0d629 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/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])) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index 5a628ed35..d4a7a7795 100644 --- a/pyemma/thermo/util/util.py +++ b/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 # ================================================================================================== diff --git a/setup.py b/setup.py index c1ee94941..3025bb7e6 100755 --- a/setup.py +++ b/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',