## Pulling the Radius of Gyration with GROMACS

We have written a little Gromacs extension that allows you
to apply a harmonic (or umbrella-like) restraint along
the radius of gyration R_{g} of a group of atoms. This
can be used for non-equilibrium pulling along R_{g}, or to
simply restrain the system along R_{g}. g_wham coming
with this code can read the tpr and pullf.xvg files, so
computing a PMF with umbrella sampling should work. Constraints along
R_{g} are not implemented. The source code is an
extension of a 4.68 development version of Gromacs,
and it is available here:

For compiling the code, please follow the instructions on the Gromacs website for compiling 4.6 versions.

PLEASE NOTE: This comes without any warranty of any
kind. We tested the code for a number of cases, and it
seems to work fine. In particular, we tested that the
width of the histograms along R_{g} during an
umbrella simulation is correct, suggesting that the force
calculations are correct. But please test the code
carefully for your purpose.

The code implements two new pull geometries, namely
`pull-geometry = radius-of-gyration` or `radius-of-gyration1`,
to be specified in the mdp file. The first, `pull-geometry
= radius-of-gyration`, indicates pulling along the normal
R_{g}, given by

R_{g}^{2} = M^{-1 }Σ_{i }m_{i} (**r**_{i} - **R**_{s})^{2},

where M = Σ_{i} m_{i} is the total
mass, and **R**_{s} = M^{-1 } Σ_{i} m_{i}**r**_{i} is the center
of mass. With `radius-of-gyration1`, the system is instead
pulled along

R_{g}^{(1)} = M^{-1 }Σ_{i} m_{i} |**r**_{i}-**R**_{s}|.

The motivation behind the R_{g}^{(1)}
coordinate is the following: When pulling along the normal
R_{g}, the forces applied to the atoms increases with the
distance of the atom from the center of mass. You see this easily when taking
the derivative of R_{g} wrt. to the atomic coordinates. Depending on
your system, this may be undesirable. When pulling along the
R_{g}^{(1)} coordinate, in contrast, the
force on the atoms does not depend on the distance from the
center of mass.

The dimensions from which the radius of gyration is
computed (x, y, and/or z) is given via the mdp option `pull-dim`. For
instance, for the R_{g} around the z-axis, use
`pull-dim = Y Y N` (meaning that only the x- and
y-coordinates are used to compute R_{g}). The normal radius of gyration
is obtained with `pull-dim = Y Y Y`.

Not surprisingly, specify the atoms used to compute
R_{g} as `pull-group1`
(or `pull-group2` etc.). Make sure to choose a good
PBC atom with `pull_pbcatom1`, such that the PBC
atom is near the center of mass of the pull
group. Otherwise, the pull code might not be able to make
the group of atoms whole, leading to a sudden jump of the
calculated R_{g} and possibly to huge pull forces.
The pull options for the force constant, pulling speed,
and initial positions work as usual. No B-states (for
free-energy perturbation) are supported.

Output: The `pullf.xvg` writes as usual the force along the
reaction coordinate. This can be used with `g_wham -if` for
computing PMFs. The `pullx.xvg` shows the center of mass
and (as the last column) the reaction coordinate. Use
`xmgrace -nxy` to see R_{g}.

We do not have a reference to cite this code. Maybe you find another paper by us that could be relevant for your work? If not, feel free to use the code, and please mention in your paper that you took the source code from our website. However, please understand that we do not maintain the code or offer other kind of support. If you find a bug, we'll be happy about a note, but please don't expect that we fix bugs on short notice. We just hope that the code is useful for someone in our community.