Analyze a 1-dimensional PMF from umbrella sampling
==================================================

| Erika McCarthy\ :sup:`1`, Şölen Ekesan\ :sup:`1`, and Darrin M.
  York\ :sup:`1`
| :sup:`1`\ Laboratory for Biomolecular Simulation Research, Institute
  for Quantitative Biomedicine and Department of Chemistry and Chemical
  Biology, Rutgers University, Piscataway, NJ 08854, USA

Learning objectives
-------------------

- Analyze umbrella sampling for a 1-dimenional reaction using NDFES to generate a PMF

Relevant literature
-------------------

- `Catalytic mechanism and pH dependence of a methyltransferase ribozyme
  (MTR1) from computational
  enzymology <https://academic.oup.com/nar/article/51/9/4508/7127216?login=true>`__
- `Surface-Accelerated String Method for Locating Minimum Free Energy
  Paths <https://pubs.acs.org/doi/full/10.1021/acs.jctc.3c01401>`__
- `Extension of the Variational Free Energy Profile and Multistate
  Bennett Acceptance Ratio Methods for High-Dimensional Potential of
  Mean Force Profile
  Analysis <https://pubs.acs.org/doi/10.1021/acs.jpca.1c00736>`__


Generating the PMF
------------------

Once the simulations are completed, take a look at the analyze_1.sh

.. code-block:: bash

      #!/bin/bash

      mkdir analysis
      mkdir analysis/dumpaves

      cd analysis

      for i in $(seq 1 32); do

        base=$(printf "img%02i" ${i})
        ndfes-PrepareAmberData.py -d ../${base}.disang -i ../${base}.dumpave -o dumpaves/${base}.dumpave -r 1 >> metafile  

      done

      ndfes_omp --mbar -w 0.15 --nboot 0 -c metafile.chk metafile

      ndfes-path_omp --chk metafile.chk --ipath metafile --npathpts 32 --nsplpts 32 --opath path


This script will create an analysis directory and use the NDFES program
to generate the potential of mean force. Each simulation generates a
dumpave file that tracks the values of the restrained properties,
including the reaction coordinate. For each dumpave file, the program
ndfes-PrepareAmberData.py will extract just the data related to the
reaction coordinate and output this to a new dumpave file in
analysis/dumpaves. It will also generate a file called metafile that
will indicate the index of the Hamiltonian (default is zero), the
temperature (default is 298), the path to the dumpave file, the umbrella
center, and the force constant.

Next, the program ndfes_omp is used to solve the MBAR equations with a
bin width equal to 0.15. Setting the --nboot flag equal to zero means
that we will perform no bootstrap resamples. If you wish to put error
bars on your data points, nboot can be set to a nonzero number like 50,
but this analysis will take longer with more resamples. The output of
the program is an xml file called metafile.chk. This contains the value
of the free energy in each sampled bin.

Finally, we use the ndfes-path_omp program to print the PMF information
to a file called path. In this example, we are merely asking it to print
the free energy along a path that is defined by the 32 control points in
the metafile (ie. where we placed the umbrella windows). To create a
smoother curve, one may increase the number of spline points to represent the path (--nsplpts). When
using interpolation, the default setting is to use radial basis
functions (rbf). Additional options will be discussed in upcoming
exercises with more complex reaction paths.

Run the analysis script and take a look at the ouput.

.. code-block:: bash

      [user@cluster] bash analyze_1.sh
      [user@cluster] ls analysis
      dumpaves  metafile  metafile.chk  path  path.rbf.0.dat  path.xml

Plot the PMF using xmgrace:

.. code-block:: bash

      [user@cluster] cd analysis
      [user@cluster] xmgrace -block path.rbf.0.dat -bxy 3:4