Merge pull request #40 from markovmodel/add_pbc

Add pbc

CHANGELOG

@@ -108,4 +108,8 @@ Version 0.2.3 (Beta):
     - Fixes: self-consistency check in cset
     - Removes several misleading warnings
 
+Version 0.2.4 (Beta):
+
+    - Allow periodicity when computing umbrella sampling bias energies
+
 Upcoming:

conda-recipe/meta.yaml

@@ -21,6 +21,9 @@ requirements:
     - msmtools >=1.1.3
 
 test:
+  source_files:
+   - test/*.py
+   - .coveragerc
   requires:
     - nose
     - coverage ==4

conda-recipe/run_test.py

@@ -5,45 +5,16 @@
 import shutil
 import re
 
-src_dir = os.getenv('SRC_DIR')
-
-def coverage_report():
-    fn = '.coverage'
-    assert os.path.exists(fn)
-    build_dir = os.getenv('TRAVIS_BUILD_DIR')
-    dest = os.path.join(build_dir, fn)
-    print( "copying coverage report to", dest)
-    shutil.copy(fn, dest)
-    assert os.path.exists(dest)
-
-    # fix paths in .coverage file
-    with open(dest, 'r') as fh:
-        data = fh.read()
-    match= '"/home/travis/miniconda/envs/_test/lib/python.+?/site-packages/.+?/(thermotools/.+?)"'
-    repl = '"%s/\\1"' % build_dir
-    data = re.sub(match, repl, data)
-    os.unlink(dest)
-    with open(dest, 'w+') as fh:
-       fh.write(data)
-
 nose_run = [
     "nosetests",
-    os.path.join(src_dir, "test"),
+    "test",
     "-vv",
     "--with-coverage",
     "--cover-inclusive",
     "--cover-package=thermotools",
     "--with-doctest",
     "--doctest-options=+NORMALIZE_WHITESPACE,+ELLIPSIS"]
 
-shutil.copyfile(
-    os.path.join(src_dir, ".coveragerc"),
-    os.path.join(os.getcwd(), ".coveragerc"))
-
 res = subprocess.call(nose_run)
 
-# move .coverage file to git clone on Travis CI
-if os.getenv('TRAVIS', False):
-    coverage_report()
-
 sys.exit(res)

ext/util/_util.c

@@ -1,7 +1,7 @@
 /*
 * This file is part of thermotools.
 *
-* Copyright 2015, 2016 Computational Molecular Biology Group, Freie Universitaet Berlin (GER)
+* Copyright 2015-2017 Computational Molecular Biology Group, Freie Universitaet Berlin (GER)
 *
 * thermotools is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
@@ -19,6 +19,12 @@
 
 #include "_util.h"
 
+inline double wrap(double x, const double width, const double inverse_width) {
+    if(width > 0.0) {
+        return x - nearbyint(x * inverse_width) * width;
+    } else return x;
+}
+
 /***************************************************************************************************
 *   sorting
 ***************************************************************************************************/
@@ -162,6 +168,7 @@ extern int _get_therm_state_break_points(int *T_x, int seq_length, int *break_po
 
 extern void _get_umbrella_bias(
     double *traj, double *umbrella_centers, double *force_constants,
+    double *width, double *inverse_width,
     int nsamples, int nthermo, int ndim, double *bias)
 {
     int i, j, s, k;
@@ -184,8 +191,8 @@ extern void _get_umbrella_bias(
                     if(0.0 == fc)
                         continue;
                     isum += fc
-                        * (traj[sdim + i] - umbrella_centers[kdim + i])
-                        * (traj[sdim + j] - umbrella_centers[kdim + j]);
+                        * wrap(traj[sdim + i] - umbrella_centers[kdim + i], width[i], inverse_width[i])
+                        * wrap(traj[sdim + j] - umbrella_centers[kdim + j], width[j], inverse_width[j]);
                 }
                 sum += isum;
             }

ext/util/_util.h

@@ -1,7 +1,7 @@
 /*
 * This file is part of thermotools.
 *
-* Copyright 2015, 2016 Computational Molecular Biology Group, Freie Universitaet Berlin (GER)
+* Copyright 2015-2017 Computational Molecular Biology Group, Freie Universitaet Berlin (GER)
 *
 * thermotools is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
@@ -77,6 +77,7 @@ extern int _get_therm_state_break_points(int *T_x, int seq_length, int *break_po
 
 extern void _get_umbrella_bias(
     double *traj, double *umbrella_centers, double *force_constants,
+    double *width, double *inverse_width,
     int nsamples, int nthermo, int ndim, double *bias);
 
 /***************************************************************************************************

ext/util/util.pyx

@@ -1,6 +1,6 @@
 # This file is part of thermotools.
 #
-# Copyright 2015, 2016 Computational Molecular Biology Group, Freie Universitaet Berlin (GER)
+# Copyright 2015-2017 Computational Molecular Biology Group, Freie Universitaet Berlin (GER)
 #
 # thermotools is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
@@ -54,6 +54,7 @@ cdef extern from "_util.h":
     # bias calculation tools
     void _get_umbrella_bias(
         double *traj, double *umbrella_centers, double *force_constants,
+        double *width, double *inverse_width,
         int nsamples, int nthermo, int ndim, double *bias)
     # transition matrix renormalization
     void _renormalize_transition_matrix(double *p, int n_conf_states, double *scratch_M)
@@ -380,7 +381,8 @@ def _overlap_post_hoc_RE(
 def get_umbrella_bias(
     _np.ndarray[double, ndim=2, mode="c"] traj not None,
     _np.ndarray[double, ndim=2, mode="c"] umbrella_centers not None,
-    _np.ndarray[double, ndim=3, mode="c"] force_constants not None):
+    _np.ndarray[double, ndim=3, mode="c"] force_constants not None,
+    _np.ndarray[double, ndim=1, mode="c"] width not None):
     r"""
     Restrict full list of samples to a subset and relabel configurational state indices.
 
@@ -402,10 +404,15 @@ def get_umbrella_bias(
     nthermo = umbrella_centers.shape[0]
     ndim = traj.shape[1]
     bias = _np.zeros(shape=(nsamples, nthermo), dtype=_np.float64)
+    inverse_width = _np.zeros(shape=(ndim,), dtype=_np.float64)
+    idx = _np.where(width > 0.0)[0]
+    inverse_width[idx] = 1.0 / width[idx]
     _get_umbrella_bias(
         <double*> _np.PyArray_DATA(traj),
         <double*> _np.PyArray_DATA(umbrella_centers),
         <double*> _np.PyArray_DATA(force_constants),
+        <double*> _np.PyArray_DATA(width),
+        <double*> _np.PyArray_DATA(inverse_width),
         nsamples,
         nthermo,
         ndim,

test/test_util.py

@@ -201,7 +201,8 @@ def test_get_umbrella_bias_binary():
     umbrella_centers[1, :] = 1.0
     force_constants = np.array([
         np.zeros(shape=(ndim, ndim), dtype=np.float64), np.eye(ndim, dtype=np.float64)])
-    bias = util.get_umbrella_bias(traj, umbrella_centers, force_constants)
+    width = np.zeros(shape=(ndim,), dtype=np.float64)
+    bias = util.get_umbrella_bias(traj, umbrella_centers, force_constants, width)
     ref = np.vstack((
         np.zeros(shape=(nsamples)),
         0.5 * ndim * np.linspace(-1.0, 1.0, nsamples)**2)).T.astype(np.float64)