Pulling the Radius of Gyration with GROMACS
We have written a 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 (typically molecule).
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} \Sigma_i m_i (\textbf{r}_i - \textbf{R}_s)^2, $$
where $M = \Sigma_i m_i$ is the total mass, and $R_s = M^{-1} \Sigma_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} \Sigma_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.
