## Pyranose ring puckering

Note:This is intentionally an advanced tutorial, showing one way of using of an approximate PMF that describes slow modes of a small molecule to setup a replica exchange simulation of that small molecule in explicit solvent. The outcome is twofold: (1) enchanced sampling of the solute degrees of freedom; and (2) accurate PMF of the slow modes involved into the setup. You will want to make reference to the following two items:

• V. Babin and C. Sagui. Conformational free energies of methyl-alpha-L-iduronic and methyl-beta-D-glucuronic acids in water. J. Chem. Phys. 132, 104018 (2010). Here is a preprint, which describes what is being computed, and why and how;
• The tar file (in gzipped form) that contains the inputs and outputs.

/////////////////////////////////////////////////////////////////////////////

The aim of this document is to present a simulation protocol that can be
used to efficiently sample puckering configurations of the
methyl-$\alpha$-L-iduronic acid in the explicit solvent bath. The
methodology is explained in this paper:

V. Babin and C. Sagui. Conformational free energies of
methyl-alpha-L-iduronic and methyl-beta-D-glucuronic acids in water.
J. Chem. Phys. 132, 104018 (2010)

The purpose of this document is to describe the technicalities. To make a
long story short: the idea is to setup a replica-exchange scheme by means
of biasing of the system along its slowest, difficult to sample, modes. In
general case the latter are not so easy to identify and capture, yet for
small compounds it is often possible to identify a suitable collective
variable'' (or a reaction coordinate'') that provides a good enough''
description of the conformational landscape. It is then possible to use
the corresponding potential of mean force (PMF), known also as Landau
free energy'', to setup an efficient replica exchange simulation. For the
molecule considered here the slow modes in question are related to the
changes of the pyranose ring puckering. These conformational changes can
be described in terms of the intraring dihedral angles O5--C1--C2--C3 and
C3--C4--C5--O5 (see Fig.2 in the preprint). This is not the only possibility,
obviously. As explained in the preprint, we proceed in stages (the step
numbers used below coincide with the folder names in the archive available here.

=============================================================================
0. Preliminaries.
-----------------------------------------------------------------------------

We assume that the reader is familiar with basic AMBER concepts and
do not repeat them here [that is, we do not explain how to create the
topology file, initial coordinates, etc].

This step is a short, regular molecular dynamics run at room temperature
in implicit solvent. The goal here is to examine the magnitude and
characteristic time scale of the thermal fluctuations of the O5--C1--C2--C3
and C3--C4--C5--O5 dihedrals. The MDIN file can be found in the
ncsu-hremd-tutorial/0/input folder of the supplementary files and
basically is nothing significantly beyond a typical GB/SA run. To get
the dihedrals we could dump the trajectory and use ptraj afterwards,
yet in our case the dihedrals are readily available through the
ncsu_abmd subsystem ran in the ANALYSIS mode. More precisely, the text
just after the &cntrl namelist instructs sander to compute the values
of the dihedrals every 100 steps and dump them into the monitor.txt''
file. For the ncsu_abmd ran in the ANALYSIS mode the latter simply
contains MD time (in picoseconds) in the first column followed the values
of the collective variables in the subsequent columns.

Once the run is complete, one can visualize the time series of the
dihedrals and conclude that the magnitude of the fluctuations is
about 0.2 radian and the pertinent time is about 0.5 ps (these values
are obviously very sensitive to the definitions; we used the old, proven
by time stare-and-pick-a-value'' approach which consists of looking at
the plots and trying to come up with some number).

=============================================================================
1. Implicit solvent.
-----------------------------------------------------------------------------

As explained in the preprint, the goal of this step is to find an
(approximate) potential of mean force (PMF, or Landau free energy) of the
O5--C1--C2--C3 and C3--C4--C5--O5 angles. We could have skipped this part
and proceeded with explicit water from the beginning, yet it is reasonable
to expect a great degree of similarity between the PMFs computed in
implicit versus explicit water models. The idea therefore is to get initial
approximation using implicit water and use it as the starting point for
more demanding explicit water computations.

We proceed in two stages: coarse "flooding", followed by finer "flooding".
The "flooding" technically amounts to changing the "mode" parameter in
the ncsu_abmd section of the MDIN file, and defining the scales: the time
scale, that controls the speed of the biasing potential buildup, and the
resolutions for the collective variables (defined separately for each of
them). There are no strong, rigorous arguments behind the choice or the
parameters. On the other hand, as we argue more in the preprint, the scales
observed in a short equilibrium run appear to give a reasonable first guess.

The values of the dihedrals, along with the value of biasing potential,
again get dumped to the "monitor.txt" file. Once all relevant minima are
visited (which in this case are: two chairs and several boats), there is
no point to continue the "coarse" simulation further and we proceed to the
"finer" run. A larger value is used for the time scale in this sub-step
and hence the dynamics is closer to equilibrium (disturbed less by the
"growing" biasing potential). The outcome is (supposedly) more accuarte
PMF (even more accurate PMF can be computed by performing an equilibrium
simulation biased by the last non-equilibrium PMF; this is not necessary
for our needs, however).

=============================================================================
2. Explicit solvent.
-----------------------------------------------------------------------------

Our ultimate goal is to have a "good enough" PMF for the dihedrals in
explicit water and use it for the replica exchange scheme. As we alluded
earlier, it is not unreasonable to expect this "explicit PMF" to be
similar to the one computed using an implicit solvent model. We therefore
use the "fine" PMF from the previous step as our starting point for
the explicit water computations.

We do "flooding" first in order to get closer to the true PMF. We use the
values of the parameters from the fine run in implicit water. As before,
we stop the simulation once everything relevant has been visited.

At that stage we hope to have a reasonable approximation of the true PMF.
In order to assess the quaility of the approximation, we do a relatively
brief (in terms of the statistics gathered) biased run at equilibrium
(umbrella subdirectory in the data folder). The approximation appears to
be good enough as the dihedrals sweep all relevant values (please see
the preprint for the more detailed presentation).

=============================================================================
3. (Hamiltonian) replica exchange in explicit solvent.
-----------------------------------------------------------------------------

Having the (approximate) PMF for the intraring dihedrals, it is easy to
setup a (Hamiltonian) replica exchange simulation that exploits the PMF
to facilitate the barriers crossing (much like the high-temperature replicas
in regular parallel tempering method do).

There are two folders holding the 'equilibration' and 'production'
files respectively. The latter is way too short as for the real
"production" (it is 3 x 1ns while something like 3 x 100ns would be
more appropriate).

We assume that the reader is familiar with doing replica-exchange
simulations using sander: basically the 'groups' file is expected to list
command line arguments to be parsed by different instances of sander;
every instance (replica) has its own MDIN file (as well as MDOUT, MDCRD,
etc). Compared to the previous step, the ncsu_abmd section is replaced
with ncsu_bbmd; this triggers a slightly different code path that does
actually perform the "exchange" part. The replica with rank 0 (first line
in the 'groups' file) has some additional parameters in the beginning of
the ncsu_bbmd section: exchange attempts frequency, name of the file to
which the exchange statistics gets dumped occasionally, and the name of
the file that stores the state of the pseudo-random generator used for
the "exchange" part (mt19937) [if the latter does not exist, the generator
can be initialized with a seed; please consult the manual for the full
reference].

-----------------------------------------------------------------------------

(c) V.Babin and C.Sagui 11/27/2009

/////////////////////////////////////////////////////////////////////////////