[gmx-developers] Restraining the radius of gyration
Albert Mao
the_one_smiley at yahoo.com
Fri Aug 11 10:20:30 CEST 2006
(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
More information about the gromacs.org_gmx-developers
mailing list