From: Aron Broom (broomsday_at_gmail.com)
Date: Thu May 24 2012 - 21:37:42 CDT
Shaon,
It does sound odd, I have one suggestion and one question:
Question: you say you don't need the *.xsc file, but you say that you take
equilibrated coordinates and then resample the velocities, what about the
periodic cell dimensions and origin? Do you just copy them from the *.xsc
file, or use VMD to recalculate them? It you just reused the values from
the start of the equilibration, they could be off.
Suggestion: Since increasing the wall to 60-70 seemed to make the wall
actually work, but then you still don't get any exploration, maybe redo
that simulation, but change "fullsamples" to something unreasonably low
like 10, that way the biasing will be applied sooner. But the second
important thing here, is that samples actually need to be collected, and if
your Rg is sitting outside your boundaries, I don't think that is
happening, which means the biasing may not be in effect (you are not
hitting fullsamples for any of the bins). You can check this by checking
the applied forces in the traj file, but you should also extend the
lowerboundary to 4, BUT then specifically specify lowerwall 7.0, that way
you are still collecting information when you are at 6.8 that will drive it
away from 6.8, but also still have the wall where you want it.
~Aron
On Thu, May 24, 2012 at 9:53 PM, Shaon Chakrabarti <shaonc_at_gmail.com> wrote:
> 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:
>
> PROBLEM 1)
> 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.
>
>
> PROBLEM 2)
> 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
> results.
> 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
> towards
> low Rg values (with or without wall potentials).
>
> Following Aron's suggestion, I started increasing the wall constants from
> the
> original value of 10 and it does indeed help--but only up to a point. At
> wall
> 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
> better...the
> 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 long
> (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
> simulations,
> 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. Using
> VMD, I select that frame and write out a pdb file. I then reinitialize the
> velocities,
> 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,
> Shaon
>
>
>
>
> On Thu, May 24, 2012 at 7:01 AM, Giacomo Fiorin <giacomo.fiorin_at_gmail.com>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 <shaonc_at_gmail.com>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
>>>
>>> VMD MEASURE RGYR:
>>> 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 <shaonc_at_gmail.com>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 <
>>>> giacomo.fiorin_at_gmail.com> 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:
>>>>>
>>>>> http://www.ks.uiuc.edu/Research/namd/2.9/ug/node55.html#SECTION0001322170000000000000
>>>>> 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 <broomsday_at_gmail.com>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 <shaonc_at_gmail.com>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
>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>
-- 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