Pulling the Radius of Gyration with GROMACS

Rg potential

We have written a little Gromacs extension that allows you to apply a harmonic (or umbrella-like) restraint along the radius of gyration Rg of a group of atoms. This can be used for non-equilibrium pulling along Rg, or to simply restrain the system along Rg. 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 Rg are not implemented. The source code is an extension of a 4.68 development version of Gromacs, and it is available here:

Download source code

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 Rg 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 Rg, given by

Rg2 = M-1 Σi mi (ri - Rs)2,

where M = Σi mi is the total mass, and Rs = M-1 Σi miri is the center of mass. With radius-of-gyration1, the system is instead pulled along

Rg(1) = M-1 Σi mi |ri-Rs|.

The motivation behind the Rg(1) coordinate is the following: When pulling along the normal Rg, 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 Rg wrt. to the atomic coordinates. Depending on your system, this may be undesirable. When pulling along the Rg(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 Rg around the z-axis, use pull-dim = Y Y N (meaning that only the x- and y-coordinates are used to compute Rg). The normal radius of gyration is obtained with pull-dim = Y Y Y.

Not surprisingly, specify the atoms used to compute Rg 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 Rg 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 Rg.

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.