[gmx-developers] Restraining the radius of gyration

[Albert Mao]

(Let me know if I should reply to just Dr. Berk Hess instead of the whole
mailing list)

Your answer was surprising to me because g_gyrate shows the radius of gyration
converging towards and oscillating around my desired value, but I can see how
the code could be incorrect yet still produce this result.

I thought that I inserted the call to restrain_rg at a point where the
coordinates x are already shifted to take periodic boundary conditions into
account.  I admit that I don't have a thorough understanding of how periodic
boundary conditions work in Gromacs.  Is there a place I can insert my call to
restrain_rg so that it will handle everything correctly, or should I modify the
function somehow, or what?

We want to use the restraint potential for a number of things, but one of them
is generating a population of starting configurations in order to calculate
transition probabilities between regions in Rg space.

Thanks for your help!
-Albert Mao

--- Berk Hess <hessb at mpip-mainz.mpg.de> wrote:
> 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.
> >  
> >
> 
> _______________________________________________
> 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.
> 


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