[gmx-developers] Restraining the radius of gyration

[Berk Hess]

The code is not correct as it does not take periodic boundary conditions 
into account.
Your group of atoms could be located on a box edge, in which case part
of the group is on one side and part of the group on the other side of 
the box.
This will result in an incorrect Rg.

I don't think this code will be useful for many people.
But what do you need it for?

Berk.

Albert Mao wrote:

>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 
>_______________________________________________
>gmx-developers mailing list
>gmx-developers at gromacs.org
>http://www.gromacs.org/mailman/listinfo/gmx-developers
>Please don't post (un)subscribe requests to the list. Use the 
>www interface or send it to gmx-developers-request at gromacs.org.
>  
>