Add unit tests for replica_exchange_EE.py, the utils module and a part of the analysis module

.circleci/config.yml

@@ -53,9 +53,13 @@ jobs:
           name: Run unit tests
           command: |
             source $HOME/pkgs/bin/GMXRC
-            pip3 install pytest 
+            pip3 install pytest
+            pip3 install pytest-mpi
             pip3 install pytest-cov
             pytest -vv --disable-pytest-warnings --cov=ensemble_md --cov-report=xml --color=yes ensemble_md/tests/
+            # COVERAGE_FILE=.coverage_1 pytest -vv --disable-pytest-warnings --cov=ensemble_md --cov-report=xml --color=yes ensemble_md/tests/
+            # COVERAGE_FILE=.coverage_2 mpirun -np 4 pytest -vv --disable-pytest-warnings --cov=ensemble_md --cov-report=xml --color=yes ensemble_md/tests/test_mpi_func.py --with-mpi
+            # coverage combine .coverage_*
 
       - run:
           name: CodeCov

.gitignore

@@ -39,6 +39,7 @@ pip-delete-this-directory.txt
 htmlcov/
 .tox/
 .coverage
+.coverage_*
 .coverage.*
 .cache
 .pytest_cache

docs/simulations.rst

@@ -222,20 +222,20 @@ include parameters for data analysis here.
       :code:`/usr/local/gromacs/bin/gmx`, the path returned by the command :code:`which gmx`) should be used.
       Note that REXEE only works with MPI-enabled GROMACS. 
 
-3.2. Input files
-----------------
+3.2. Input settings
+-------------------
 
   - :code:`gro`: (Required)
-      The input system configuration in the form of GRO file(s) used to initiate the REXEE simulation. If only one GRO file is specified,
+      The path of the input system configuration in the form of GRO file(s) used to initiate the REXEE simulation. If only one GRO file is specified,
       it will be used to initiate all the replicas. If multiple GRO files are specified (using the YAML syntax),
       the number of GRO files has to be the same as the number of replicas. 
   - :code:`top`: (Required)
-      The input system topology in the form of TOP file(s) used to initiate the REXEE simulation. If only one TOP file is specified,
+      The path of the input system topology in the form of TOP file(s) used to initiate the REXEE simulation. If only one TOP file is specified,
       it will be used to initiate all the replicas. If multiple TOP files are specified (using the YAML syntax),
       the number of TOP files has to be the same as the number of replicas. In the case where multiple TOP and GRO files are specified,
       the i-th TOP file corresponds to the i-th GRO file.
   - :code:`mdp`: (Required)
-      The input MDP file used to initiate the REXEE simulation. Specifically, this input MDP file will serve as a template for
+      The path of the input MDP file used to initiate the REXEE simulation. Specifically, this input MDP file will serve as a template for
       customizing MDP files for all replicas. Therefore, the MDP template must have the whole range of :math:`λ` values. 
       and the corresponding weights (in fixed-weight simulations). This holds for REXEE simulations for multiple serial mutations as well.
       For example, in an REXEE simulation that mutates methane to ethane in one replica and ethane to propane in the other replica, if

ensemble_md/analysis/analyze_matrix.py

@@ -16,6 +16,7 @@
 from matplotlib import cm
 from ensemble_md.utils.exceptions import ParseError
 from ensemble_md.utils.exceptions import ParameterError
+from ensemble_md.analysis import synthesize_data
 
 
 def calc_transmtx(log_file, expanded_ensemble=True):
@@ -122,23 +123,30 @@ def calc_equil_prob(trans_mtx):
     return equil_prob
 
 
-def calc_spectral_gap(trans_mtx, atol=1e-8):
+def calc_spectral_gap(trans_mtx, atol=1e-8, n_bootstrap=50, bootstrap_seed=None):
     """
-    Calculates the spectral gap of the input transition matrix.
+    Calculates the spectral gap of the input transition matrix and estimates its
+    uncertainty using the bootstrap method.
 
     Parameters
     ----------
     trans_mtx : np.ndarray
-        The input state transition matrix
+        The input transition matrix
     atol: float
         The absolute tolerance for checking the sum of columns and rows.
+    n_bootstrap: int
+        The number of bootstrap iterations for uncertainty estimation.
+    bootstrap_seed: int
+        The seed for the random number generator for the bootstrap method.
 
     Returns
     -------
     spectral_gap : float
-        The spectral gap of the input transitio n matrix
+        The spectral gap of the input transitio n matrix.
+    spectral_gap_err : float
+        The estimated uncertainty of the spectral gap.
     eig_vals : list
-        The list of eigenvalues
+        The list of eigenvalues.
     """
     check_row = sum([np.isclose(np.sum(trans_mtx[i]), 1, atol=atol) for i in range(len(trans_mtx))])
     check_col = sum([np.isclose(np.sum(trans_mtx[:, i]), 1, atol=atol) for i in range(len(trans_mtx))])
@@ -159,7 +167,58 @@ def calc_spectral_gap(trans_mtx, atol=1e-8):
 
     spectral_gap = np.abs(eig_vals[0]) - np.abs(eig_vals[1])
 
-    return spectral_gap, eig_vals
+    # Estimate the uncertainty of the spectral gap
+    spectral_gap_list = []
+    n_performed = 0
+    while n_performed < n_bootstrap:
+        mtx_boot = synthesize_data.synthesize_transmtx(trans_mtx, seed=bootstrap_seed)[0]
+        check_row_boot = sum([np.isclose(np.sum(mtx_boot[i]), 1, atol=atol) for i in range(len(mtx_boot))])
+        check_col_boot = sum([np.isclose(np.sum(mtx_boot[:, i]), 1, atol=atol) for i in range(len(mtx_boot))])
+        if check_row_boot == len(mtx_boot):
+            eig_vals_boot, _ = np.linalg.eig(mtx_boot.T)
+        elif check_col_boot == len(mtx_boot):
+            eig_vals_boot, _ = np.linalg.eig(mtx_boot)
+        else:
+            continue
+
+        n_performed += 1
+        eig_vals_boot = np.sort(eig_vals_boot)[::-1]
+        spectral_gap_list.append(np.abs(eig_vals_boot[0]) - np.abs(eig_vals_boot[1]))
+
+    spectral_gap_err = np.std(spectral_gap_list, ddof=1)
+
+    return spectral_gap, spectral_gap_err, eig_vals
+
+
+def calc_t_relax(spectral_gap, exchange_period, spectral_gap_err=None):
+    """
+    Calculates the relaxation time given the spectral gap of a transition matrix of interest.
+    By defintion, the relaxation time is equal to the exchange period divided by the spectral gap.
+
+    Parameters
+    ----------
+    spectral_gap: float
+        The input spectral gap.
+    exchange_period : float
+        The exchange period of the simulation in ps.
+    spectral_gap_err : float
+        The uncertainty of the spectral gap, which is used to calculate the uncertainty of the relaxation time using
+        error propagation.
+
+    Returns
+    -------
+    t_relax : float
+        The relaxation time in ps.
+    t_relax_err : float
+        The uncertainty of the relaxation time in ps.
+    """
+    t_relax = exchange_period / spectral_gap
+    t_relax_err = None
+
+    if spectral_gap_err is not None:
+        t_relax_err = exchange_period * spectral_gap_err / spectral_gap ** 2  # error propagation
+
+    return t_relax, t_relax_err
 
 
 def split_transmtx(trans_mtx, n_sim, n_sub):