[gmx-developers] Restraining the radius of gyration

[Albert Mao]

Greetings!

I have implemented the ability to impose a restraining potential on the radius
of gyration of a group of atoms in Gromacs 3.3.1.  The definition of "radius of
gyration" I'm using here is:
Rg = sqrt(SumOverParticles(Mass_i * |X_i - X_cm|^2) / TotalMass)
The potential is:
V = 1/2 k (Rg - TargetRg)^2
The user specifies the parameters of the constraint in the .mdp file:
user1-grps specifies the group(s) of atoms to be restrained
userreal1 specifies TargetRg, the desired radius of gyration (nanometers)
userreal2 specifies k, the force constant (kJ / mol / nm^2)

I am posting my (very simple) code here in order to ask two questions:
- Would many other people find this capability useful?  If so, somebody with
CVS privileges may want to integrate my code into the main codebase, but only
after some additional engineering (for example, the .mdp and .tpr formats
should be modified instead of using the user-defined parameters).
- Does the code look correct?  I'm pretty confident about the math, but I had
to familiarize myself with the Gromacs source code in a short amount of time
and could have missed something (for example, I'm pretty sure the center of
mass gets calculated elsewhere, but I just recalculate it within my function
for simplicity).

Thanks!
-Albert Mao
Washington University in St. Louis
Medical Scientist Training Program

---------- Added to mdlib/sim_util.c near the end of do_force, immediately
before the call to calc_f_el ----------
    /* Computes the radius of gyration restraints. */
    if (inputrec->userreal1 > 0 || inputrec->userreal2 > 0)
        restrain_rg(NULL, start, homenr, mdatoms, x, f,
        inputrec->userreal1, inputrec->userreal2);

---------- Added to mdlib/sim_util.c ----------
static void restrain_rg(FILE *fp, int start, int homenr, t_mdatoms *mdatoms,
                        rvec x[], rvec f[], real rg, real k)
{
    int i, d;
    real totalMass = 0;
    real bigI = 0;
    rvec cm; for (d = 0; d < DIM; ++d) cm[d] = 0;
    real mx;
    real deviation;
    real forceFactor;
    
    /* For now, this loop redundantly calculates totalMass, cm, and bigI on all
nodes.
     * It could be parallelized by summing over only the atoms belonging to a
node
     * and accumulating the sums around the ring.
     */
    for (i = 0; i < mdatoms->nr; ++i) if (mdatoms->cU1[i] == 0) {
        totalMass += mdatoms->massA[i];
        for (d = 0; d < DIM; ++d) {
            mx = mdatoms->massA[i] * x[i][d];
            cm[d] += mx;
            bigI += mx * x[i][d];
        }
    }
    
    if (totalMass > 0) {
        for (d = 0; d < DIM; ++d) {
            cm[d] /= totalMass;
            bigI -= totalMass * cm[d] * cm[d];
        }
        deviation = sqrt(bigI/totalMass) - rg;
        forceFactor = -k * deviation / sqrt(bigI*totalMass);
        for (i = start; i < start+homenr; ++i) if (mdatoms->cU1[i] == 0) {
            for (d = 0; d < DIM; ++d)
                f[i][d] += forceFactor * mdatoms->massA[i] * (x[i][d] - cm[d]);
        }
    }
}

__________________________________________________
Do You Yahoo!?
Tired of spam?  Yahoo! Mail has the best spam protection around 
http://mail.yahoo.com