Re: Issues with 'Gyration' in Colvars module

From: Shaon Chakrabarti (
Date: Thu May 24 2012 - 20:53:25 CDT

Hi Aron, Giacomo
Thanks so much for all the suggestions...Problem 1) is almost
resolved, but Problem 2) is still creating problems.
Here's what I find:

I defined Rg using the mass weighting definition as follows:

Rg = sqrt[ 1/M * sum over all N atoms
{m*(x-COMx)^2+m*(y-COMy)^2+m*(z-COMz)^2} ]
and similarly for the COG case.

Now, Rg (COM) is 5.5785031 while Rg (COG) is 5.5799608

So indeed, Rg (COG) defined with the mass-weighted formula matches
the colvars.traj result, as Giacomo had said. However, the formula
most closely matching the VMD rgyr result seems to be the one
with COG (without mass-weighting)! I will check once again to see
if I'm making a mistake however.

This is still causing problems. First let me mention that my system
is a 25-mer alkane in water, so though it does have a global minimum
at around Rg = 4 angstroms, it is expected to collapse from an extended
state on a timescale of 1 ns. Without any biases, that's precisely what I
observe in my simulation. I've applied HARMONIC biases (with spring
constant as low as 1 kcal/mol.A^2) to the Rg colvar with perfectly fine
I have also used ABF with the 'distance' component (order parameter
end-to-end distance) and the results are fine, even with wall constant 10.
The moment I use ABF along with Rg, the system tends to get rapidly pushed
low Rg values (with or without wall potentials).

Following Aron's suggestion, I started increasing the wall constants from
original value of 10 and it does indeed help--but only up to a point. At
constants of about 60-70, I observe that instead of hovering around 4
angstroms(like before), the system now stays at around 6.8-7 angstroms.
(Remember, I had my lower and upper boundaries at 7 and 11). Unfortunately,
even after increasing the wall constant to 5x10^5, this doesn't get much
system fluctuates in the region of 7.0-7.4 angstroms and does not explore
higher values of Rg even in a 100000 step (0.2 ns) run. I will try to run a
(3-4 ns) trajectory to make sure the system does indeed stay around 7 all
the time. If I increase the wall constants even more, I start getting
RATTLE errors.

I am using constant P,T simulations, but since I'm not restarting the
I don't need the *.xsc file. From my original long equilibrated trajectory
(with no
biases), I choose a configuration where Rg lies between 7 and 11 angstroms.
VMD, I select that frame and write out a pdb file. I then reinitialize the
apply the necessary bias and generate the trajectory.

So at this point it seems that ONLY when ABF is used with Rg, the system
behaves weird.

Thanks again for the help,

On Thu, May 24, 2012 at 7:01 AM, Giacomo Fiorin <>wrote:

> Hello Shaon, I don't see the masses in your formula, what happens if you
> include them? I have a feeling that the numbers would coincide in that
> case.
> On the other question, Aron wrote you already what would have been my
> thoughts on the trajectory.
> Giacomo
> On Wed, May 23, 2012 at 6:25 PM, Shaon Chakrabarti <>wrote:
>> ....And here are the details of 1)
>> The following calculations have been
>> done on a pdb file (which I could send
>> separately) grabbed from the end of
>> a trajectory run. The psf file I used in
>> VMD is the same one I used to run
>> the NAMD simulation.
>> MY CALCULATION: (x,y,z denote atom coordinates)
>> a) Rg defined with center of mass:
>> Rg = sqrt[1/N * sum over all N atoms{(x-COMx)^2+(y-COMy)^2+(z-COMz)^2} ]
>> With this definition I get Rg = 5.7555
>> b) Rg defined with center of geometry:
>> COGx = (1/N * sum of all x components). Similarly for COGy, COGz
>> Rg = sqrt[1/N * sum over all N
>> atoms{(x-COGx)^2+(y-COGy)^2+(z-COGz)^2} ]
>> With this definition I get Rg = 5.754163
>> this gives Rg = 5.7541623 --almost identical to the result of part b)
>> above
>> COLVARS OUTPUT (from colvars.traj file):
>> this gives Rg = 5.579920
>> Hope this helps...
>> Thanks,
>> Shaon
>> On Wed, May 23, 2012 at 3:11 PM, Shaon Chakrabarti <>wrote:
>>> Hi Aron and Giacomo,
>>> Thanks a lot for such quick replies! Actually I had
>>> already done what Aron suggested: I first ran a 6ns
>>> long equilibration trajectory without the ABF module,
>>> but with gyration on.
>>> From the last parts of this trajectory, I created pdb
>>> and psf files of the molecule with Rg within the
>>> required range. Then I reinitialized the velocities and
>>> ran the trajectory with ABF included. I'm attaching the
>>> initial bit of colvars.traj file here--the range is 7-11 (i had the
>>> file for this range). The starting Rg value is about 8...
>>> so that is between 7 and 11...yet, the Rg value doesn't
>>> stay within the required range!
>>> I'm also attaching the colvars configuration file if that helps.
>>> About 1) ...I will provide the numbers and details of my code in
>>> a little while.
>>> Thanks again for all the help,
>>> Shaon
>>> On Wed, May 23, 2012 at 1:27 PM, Giacomo Fiorin <
>>>> wrote:
>>>> Hi Shaon, Aron has already answered well on 2).
>>>> About 1), can you provide the exact numbers, and more details regarding
>>>> your own calculation? Do you load the same PSF topology and coordinates in
>>>> VMD or use the masses guessed by VMD?
>>>> In the colvars module, gyration is computed as the mass-weighted square
>>>> distance from the center of geometry, taken under square root. But I
>>>> actually realize now that the formula in the user's guide:
>>>> does not include the weighting masses. Would that explain the small
>>>> difference that you saw? If so, I apologize for the omission of the masses
>>>> in the user's guide.
>>>> The gyration defined in VMD is also mass-weighted, but uses the center
>>>> of mass.
>>>> In gyration, it seems contradictory to have mass-weighted gyration
>>>> around the center of geometry and not mass, but we decided towards that
>>>> considered that gyration can be and is often used in combination with RMSD.
>>>> On Wed, May 23, 2012 at 3:56 PM, Aron Broom <>wrote:
>>>>> Hi Shaon,
>>>>> I've no idea about 1), but for 2), what happens if you specifiy Rg to
>>>>> be within 4-7 rather than 5-7? From my understanding the ABF method needs
>>>>> to sample forces before it starts to apply a bias, so if you're system
>>>>> starts outside the biasing range, it will never collect any samples, and
>>>>> hence never start to apply a bias. That would be my first thought anyway.
>>>>> If you really want to do 5-7 angstroms, you should apply a harmonic
>>>>> potential on Rg and sit it ~6 angstroms in order to equilibrate the system
>>>>> before running ABF.
>>>>> ~Aron
>>>>> On Wed, May 23, 2012 at 3:48 PM, Shaon Chakrabarti <>wrote:
>>>>>> Hi NAMD experts,
>>>>>> I'm facing two problems with the gyration component in
>>>>>> the Colvars module.
>>>>>> 1) The NAMD user's guide says that the radius of gyration
>>>>>> is defined using the Center of Geometry. However, when
>>>>>> I write my own code (using the COG definition) to compute
>>>>>> the Rg of a molecule, the Rg value is always about 0.2
>>>>>> angstroms more compared to the output in the 'colvars.traj'
>>>>>> file. The values from my code exactly match the values
>>>>>> I get if I use the 'measure rgyr' command in VMD. I also
>>>>>> computed Rg using Center of Mass, but that makes the
>>>>>> difference even larger. I'd be grateful if someone can tell
>>>>>> me if I'm making a mistake somewhere, or if Rg is
>>>>>> actually defined differently in the Colvars module.
>>>>>> P.S. My molecule (about 30 angstroms long when fully
>>>>>> extended) is in the middle of a water box of side 55
>>>>>> angstroms.
>>>>>> 2) The second problem is similar (as I understand from
>>>>>> the threads in 2008) to what Subramanian
>>>>>> Vaitheeswaran faced while using the gyration
>>>>>> component with ABF. If I specify the Rg to be (say)
>>>>>> within say 5 and 7 angstroms, and use ABF to get
>>>>>> the PMF, the system hovers at Rg value around
>>>>>> 4.5-4.8 most of the time! A final resolution doesn't
>>>>>> seem to have been posted back in 2008, and it
>>>>>> would be of immense help if someone could tell
>>>>>> me what's going on....if the problem has been
>>>>>> resolved. I tried switching the hideJacobian option
>>>>>> to yes, and increasing the lower and upper wall
>>>>>> constants, but these don't make any difference.
>>>>>> It would be great if someone could help with these two issues.
>>>>>> Thanks,
>>>>>> Shaon
>>>>> --
>>>>> Aron Broom M.Sc
>>>>> PhD Student
>>>>> Department of Chemistry
>>>>> University of Waterloo

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:21:33 CST