Implement the workflow pipeline

.gitignore

@@ -2,3 +2,6 @@
 .vscode
 *.DS_Store
 build
+.idea/
+docs/_build/
+src/alchemlyb.egg-info/

docs/index.rst

@@ -57,6 +57,7 @@ Contributions are very welcome. If you have bug reports or feature requests or q
     :caption: User Documentation
 
     install
+    workflow
     parsing
     preprocessing
     estimators

docs/visualisation.rst

@@ -12,6 +12,7 @@ visualisation tools to help user to judge the estimate.
     plot_mbar_overlap_matrix
     plot_ti_dhdl
     plot_dF_state
+    plot_convergence
 
 .. _plot_overlap_matrix:
 
@@ -131,6 +132,53 @@ Will give a plot looks like this
    A bar plot of the free energy differences evaluated between pairs of adjacent
    states via several methods, with corresponding error estimates for each method.
 
+.. _plot_convergence:
+
+Forward and Backward Convergence
+--------------------------------
+One way of determining the simulation end point is to plot the forward and
+backward convergence of the estimate using
+:func:`~alchemlyb.visualisation.plot_convergence`.
+
+Note that this is just a plotting function to plot [Klimovich2015]_ style
+convergence plot. The user need to provide the forward and backward data list
+and the corresponding error. ::
+
+    >>> import pandas as pd
+    >>> from alchemtest.gmx import load_benzene
+    >>> from alchemlyb.parsing.gmx import extract_u_nk
+    >>> from alchemlyb.estimators import MBAR
+
+    >>> bz = load_benzene().data
+    >>> data_list = [extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']]
+    >>> forward = []
+    >>> forward_error = []
+    >>> backward = []
+    >>> backward_error = []
+    >>> num_points = 10
+    >>> for i in range(1, num_points+1):
+    >>>     # Do the forward
+    >>>     slice = int(len(data)/num_points*i)
+    >>>     u_nk_coul = pd.concat([data[:slice] for data in data_list])
+    >>>     estimate = MBAR().fit(u_nk_coul)
+    >>>     forward.append(estimate.delta_f_.iloc[0,-1])
+    >>>     forward_error.append(estimate.d_delta_f_.iloc[0,-1])
+    >>>     # Do the backward
+    >>>     u_nk_coul = pd.concat([data[-slice:] for data in data_list])
+    >>>     estimate = MBAR().fit(u_nk_coul)
+    >>>     backward.append(estimate.delta_f_.iloc[0,-1])
+    >>>     backward_error.append(estimate.d_delta_f_.iloc[0,-1])
+
+    >>> from alchemlyb.visualisation import plot_convergence
+    >>> ax = plot_convergence(forward, forward_error, backward, backward_error)
+    >>> ax.figure.savefig('dF_t.pdf')
+
+Will give a plot looks like this
+
+.. figure:: images/dF_t.png
+
+   A convergence plot of showing that the forward and backward has converged
+   fully.
 
 .. [Klimovich2015] Klimovich, P.V., Shirts, M.R. & Mobley, D.L. Guidelines for
    the analysis of free energy calculations. J Comput Aided Mol Des 29, 397–411

docs/visualisation/alchemlyb.visualisation.plot_convergence.rst

@@ -0,0 +1,19 @@
+.. _visualisation_plot_convergence:
+
+Plot the Forward and Backward Convergence
+=========================================
+
+The function :func:`~alchemlyb.visualisation.plot_convergence` allows
+the user to visualise the convergence by plotting the free energy change
+computed using the equilibrated snapshots between the proper target time frames
+in both forward (data points are stored in `forward` and `forward_error`) and
+reverse (data points are stored in `backward` and `backward_error`) directions.
+The unit in the y axis could be labelled to other units by setting *units*,
+which by default is kBT. The user can pass :class:`matplotlib.axes.Axes` into
+the function to have the convergence drawn on a specific axes.
+
+Please check :ref:`How to plot convergence <plot_convergence>` for usage.
+
+API Reference
+-------------
+.. autofunction:: alchemlyb.visualisation.plot_convergence

docs/visualisation/alchemlyb.visualisation.plot_dF_state.rst

@@ -25,7 +25,7 @@ The figure could be plotted in *portrait* or *landscape* mode by setting the
 The user could pass a list of strings to `labels` to name the
 :class:`~alchemlyb.estimators` or a list of strings to `colors` to color
 the estimators differently. The unit in the y axis could be labelled to other
-units by setting `units`, which by default is kcal/mol.
+units by setting `units`, which by default is kBT.
 
 Please check :ref:`How to plot dF states <plot_dF_states>` for a complete
 example.

docs/visualisation/alchemlyb.visualisation.plot_ti_dhdl.rst

@@ -11,7 +11,7 @@ When custom labels are desirable, the user could pass a list of strings to the
 *labels* for labelling each alchemical transformation differently. The color of
 each alchemical transformation could also be set by passing a list of color
 string to the *colors*. The unit in the y axis could be labelled to other units
-by setting *units*, which by default is kcal/mol. The user can pass
+by setting *units*, which by default is kBT. The user can pass
 :class:`matplotlib.axes.Axes` into the function to have the dhdl drawn on a
 specific axes.
 

docs/workflow.rst

@@ -0,0 +1,109 @@
+Automatic workflow
+==================
+Though **alchemlyb** is a library offering great flexibility in deriving free
+energy estimate, it also provide a easy pipeline that is similar to
+`Alchemical Analysis <https://github.com/MobleyLab/alchemical-analysis>`_ and a
+step-by-step version that allows more flexibility.
+
+Fully Automatic analysis
+------------------------
+A interface similar to
+`Alchemical Analysis <https://github.com/MobleyLab/alchemical-analysis>`_
+could be excuted with a single line of command. ::
+
+    >>> import os
+    >>> from alchemtest.gmx import load_ABFE
+    >>> from alchemlyb.workflows import ABFE
+    >>> # Obtain the path of the data
+    >>> dir = os.path.dirname(load_ABFE()['data']['complex'][0])
+    >>> print(dir)
+    'alchemtest/gmx/ABFE/complex'
+    >>> workflow = ABFE(units='kcal/mol', software='Gromacs', dir=dir,
+    >>>                 prefix='dhdl', suffix='xvg', T=298, skiptime=10,
+    >>>                 uncorr='dhdl', threshold=50,
+    >>>                 methods=('mbar', 'bar', 'ti'), out='./',
+    >>>                 resultfilename='result.out', overlap='O_MBAR.pdf',
+    >>>                 breakdown=True, forwrev=10, log='result.log')
+
+This would give the free energy estimate using all of
+:class:`~alchemlyb.estimators.TI`, :class:`~alchemlyb.estimators.BAR`,
+:class:`~alchemlyb.estimators.MBAR` and the result will be written to the text
+file `result.out`. ::
+
+    ------------ --------------------- --------------------- ---------------------
+       States         MBAR (kcal/mol)        BAR (kcal/mol)         TI (kcal/mol)
+    ------------ --------------------- --------------------- ---------------------
+       0 -- 1         0.041  +-  0.001      0.041  +-  0.001      0.041  +-  0.001
+       1 -- 2         0.056  +-  0.001      0.055  +-  0.001      0.056  +-  0.001
+       2 -- 3         0.082  +-  0.001      0.082  +-  0.002      0.083  +-  0.002
+    ...
+      26 -- 27        0.766  +-  0.007      0.768  +-  0.010      0.770  +-  0.010
+      27 -- 28        0.694  +-  0.008      0.691  +-  0.011      0.690  +-  0.010
+      28 -- 29        0.620  +-  0.010      0.616  +-  0.011      0.625  +-  0.011
+    ------------ --------------------- --------------------- ---------------------
+         coul:        6.290  +-  0.021      6.168  +-  0.026      6.168  +-  0.030
+          vdw:       13.872  +-  0.061     13.852  +-  0.037     13.877  +-  0.066
+       bonded:        1.469  +-  0.009      1.447  +-  0.003      1.461  +-  0.013
+        TOTAL:       21.631  +-  0.064     21.467  +-  0.054     21.506  +-  0.074
+
+The :ref:`overlay matrix for the MBAR estimator <plot_overlap_matrix>` will be
+plotted and saved to `O_MBAR.pdf`.
+
+The :ref:`dHdl for TI <plot_TI_dhdl>` will be plotted to `dhdl_TI.pdf`.
+
+The :ref:`dF states <plot_dF_states>` will be plotted to `dF_state.pdf` in
+portrait model and `dF_state_long.pdf` in landscape model.
+
+The forward and backward convergence will be plotted to `dF_t.pdf` using
+:class:`~alchemlyb.estimators.MBAR`.
+
+.. currentmodule:: alchemlyb.workflows
+
+.. autoclass:: ABFE
+   :noindex:
+
+Semi-automatic analysis
+-----------------------
+The same analysis could also performed in steps allowing access and modification
+to the data generated at each stage of the analysis. ::
+
+    >>> import os
+    >>> from alchemtest.gmx import load_ABFE
+    >>> from alchemlyb.workflows import ABFE
+    >>> # Obtain the path of the data
+    >>> dir = os.path.dirname(load_ABFE()['data']['complex'][0])
+    >>> print(dir)
+    'alchemtest/gmx/ABFE/complex'
+    >>> # Load the data
+    >>> workflow = ABFE(software='Gromacs', dir=dir,
+    >>>                 prefix='dhdl', suffix='xvg', T=298, out='./',
+    >>>                 log='result.log')
+    >>> # Set the unit.
+    >>> workflow.update_units('kcal/mol')
+    >>> # Decorrelate the data.
+    >>> workflow.preprocess(skiptime=10, uncorr='dhdl', threshold=50)
+    >>> # Run the estimator
+    >>> workflow.estimate(methods=('mbar', 'bar', 'ti'))
+    >>> # write the result
+    >>> workflow.write(resultfilename='result.out')
+    >>> # Plot the overlap matrix
+    >>> workflow.plot_overlap_matrix(overlap='O_MBAR.pdf')
+    >>> # Plot the dHdl for TI
+    >>> workflow.plot_ti_dhdl(dhdl_TI='dhdl_TI.pdf')
+    >>> # Plot the dF states
+    >>> workflow.plot_dF_state(dF_state='dF_state.pdf')
+    >>> # Convergence analysis
+    >>> workflow.check_convergence(10, dF_t='dF_t.pdf')
+
+
+.. currentmodule:: alchemlyb.workflows.ABFE
+
+.. autofunction:: update_units
+.. autofunction:: preprocess
+.. autofunction:: estimate
+.. autofunction:: write
+.. autofunction:: plot_overlap_matrix
+.. autofunction:: plot_ti_dhdl
+.. autofunction:: plot_dF_state
+.. autofunction:: check_convergence
+

src/alchemlyb/__init__.py

@@ -1,4 +1,3 @@
-
 from ._version import get_versions
 __version__ = get_versions()['version']
 del get_versions

src/alchemlyb/constants.py

@@ -0,0 +1,5 @@
+'''Physical and mathematical constants and units.'''
+# Taken from scipy.constants since py2 doesn't support it
+k = 1.380649e-23
+N_A = 6.02214076e+23
+kJ2kcal = 0.239006

src/alchemlyb/estimators/mbar_.py

@@ -48,7 +48,9 @@ class MBAR(BaseEstimator):
     """
 
     def __init__(self, maximum_iterations=10000, relative_tolerance=1.0e-7,
-                 initial_f_k=None, method='hybr', verbose=False):
+                 initial_f_k=None, method='adaptive', verbose=False):
+        # method='adaptive' is used as it is more stable
+        # https://github.com/choderalab/pymbar/issues/419
 
         self.maximum_iterations = maximum_iterations
         self.relative_tolerance = relative_tolerance

src/alchemlyb/estimators/ti_.py

@@ -105,3 +105,46 @@ def fit(self, dHdl):
         self.states_ = means.index.values.tolist()
 
         return self
+
+    def separate_dhdl(self):
+        """
+        For transitions with multiple lambda, this function will separate the
+        dhdl with multiple columns into a list of Dataframe with a single column
+        (single lambda).
+
+        Returns
+        ----------
+        dHdl_list : list
+            A list of Series such that dHdl_list[k][n] is the potential
+            energy gradient with respect to lambda for each configuration n and
+            lambda k.
+
+        """
+        if len(self.dhdl.index.names) == 1:
+            # If only one column is present convert to series
+            assert len(self.dhdl.columns) == 1
+            name = self.dhdl.columns[0]
+            return [self.dhdl[name], ]
+        else:
+            dhdl_list = []
+            # get the lambda names
+            l_types = self.dhdl.index.names
+            # obtain bool of changed lambdas between each state
+            lambdas = self.dhdl.reset_index()[l_types]
+            diff = lambdas.diff().to_numpy(dtype='bool')
+            # diff will give the first row as NaN so need to fix that
+            diff[0, :] = diff[1, :]
+            # Make sure that the start point is set to true as well
+            diff[:-1, :] = diff[:-1, :] | diff[1:, :]
+            for i in range(len(l_types)):
+                if any(diff[:,i]) == False:
+                    # Skip if not pertubed
+                    pass
+                else:
+                    new = self.dhdl.iloc[diff[:,i], i]
+                    # drop all other index
+                    for l in l_types:
+                        if l != l_types[i]:
+                            new = new.reset_index(l, drop=True)
+                    dhdl_list.append(new)
+            return dhdl_list

src/alchemlyb/preprocessing/subsampling.py

@@ -14,7 +14,6 @@ def _check_multiple_times(df):
     else:
         return df.sort_index(0).reset_index('time').duplicated('time').any()
 
-
 def _check_sorted(df):
     return df.reset_index(0)['time'].is_monotonic_increasing
 

src/alchemlyb/tests/test_ti_estimators.py

@@ -136,3 +136,8 @@ class TestTI(TIestimatorMixin):
     def X_delta_f(self, request):
         get_dHdl, E, dE = request.param
         return get_dHdl(), E, dE
+
+def test_TI_separate_dhdl():
+    dHdl = gomc_benzene_dHdl()
+    estimator = TI().fit(dHdl)
+    assert all([isinstance(dhdl, pd.Series) for dhdl in estimator.separate_dhdl()])