From: Giacomo Fiorin (gfiorin_at_seas.upenn.edu)
Date: Mon Jul 27 2009 - 11:52:44 CDT
Hi Ali, your message was quite long and detailed, and took some time
to read. Perhaps it is better to ask more specific questions on this
mailing list, and separate them if needed. Some people may be able to
answer the questions at the end of the message, but won't get there if
what comes before is not within their current experience.
That said, I think that the behavior you got pretty normal: some
collective variables work, some others don't.
I'm glad you succeeded with the RMSD colvar: please, if you use the
obtained PMF as an estimate for the folding free energy, make sure
that you have sampled enough different configurations for every point
of the PMF (i.e. that data are uncorrelated). 10 ns might be enough,
but very frequently they are quite short given the relaxation times
involved for this type of degrees of freedom (fractions of ns).
System forces always fluctuate a lot: it's their *average* that
matters, and has to be matched by the applied forces. If you run the
same example with ABF, you'll know immediately what I'm talking about.
Bests
Giacomo
---- ----
Giacomo Fiorin
Center for Molecular Modeling at
University of Pennsylvania
231 S 34th Street, Philadelphia, PA 19104-6323
phone: (+1)-215-573-4773
fax: (+1)-215-573-6233
mobile: (+1)-267-324-7676
mail: giacomo.fiorin_<at>_gmail.com
web: http://www.cmm.upenn.edu/
---- ----
On Wed, Jul 22, 2009 at 10:18 PM, Ali Emileh<ali.emileh_at_gmail.com> wrote:
> Dear Giacomo,
>
> Thanks to your help, I finally folded the helix. Although
> I didn't use the same colvar (alpha).
>
> I tried using the alpha colvar with a couple of different HillWeight
> (changing stepwise, almost an order of magnitude from 0.04 to 0.5)
> but wasn't able to fold the helix within the 10 ns time frame
> I was hoping for. So next I tried your other suggestion to
> use harmonic steering as a works/not works solution. I
> tied the same atoms to the harmonic restraint with a
> forceConstant of 10 (colvar width was set to 0.05) and
> put a lowerWallConstant of 200 on
> the lowerBoundary (0.3; the colvar starts at 0.445). Center
> of the restraint was set to 0.5 and I asked NAMD to drive
> the colvar to 1.0. During the first 5000 steps, the colvar
> goes down and past 0.3 and then stays below that (due
> to the wall? but why did it pass the wall in the 1st place?) for
> the whole duration of the simulation and somehow
> stabilizes around 0.2. Looking at the DCD, I can somewhat
> see the partially folded helix being a bit unfolded after
> the 10 ns run. Here's my config file:
>
> =====
> colvarsTrajFrequency 1000
>
> colvar {
> name a1
> width 0.05
> outputVelocity on
> outputAppliedForce on
> #outputSystemForce on
> lowerBoundary 0.3
> upperBoundary 1.0
> lowerWallConstant 200
>
> alpha {
> residueRange 100-116
> psfSegID G
> }
> }
>
> harmonic {
> colvars a1
> forceConstant 10
> centers 0.5
> targets 1.0
> targetsNumSteps 2450000
> }
> =====
>
> Now I have been able to fold this excised helix using
> TMD and in the context of its parent protein before,
> albeit using much much larger force constants, within
>>5 ns and don't know why it's not working here. Maybe
> I'm not using enough force? Also why is that the peptide
> passes the lowerWall? I'll try to repeat this run with a
> larger force constant some time later.
>
> As with metadynamics :
> I used metadynamics with hillWieght of 0.5 and
> newHillFrequencies of 100 (default) on different colvars:
> * 2 independent dihedral colvars, one representing all the
> psi angles of the backbone and the other, all the phi
> angels (obviously with a lot of overlap in the two colvars),
> no walls.
> * 2 independent colvars of hBond and distance, again
> between all the target atoms (e.g. hbonds of NH1.109-O.105
> and distance between the same atoms). no walls.
>
> No results from the above two colvars.
>
> My first -unintelligent!- guess would be to combine the
> independent colvars linearly into one large one and apply metadynamics
> to that, but is it possible to do so with colvars of different type
> and what would be the boundaries? Plus, I can't
> physically imagine the end result of such a combination (even
> a multi-dihedral colvar) when you look at the PMF vs. colvar
> graph.
>
> Next I tried the RMSD colvar. To my delight, this one worked
> and successfully folded the helix. My first shot with 0.04
> hillWeight and newHillFrequencies of 10 failed but I reduced
> the frequency to 100 (default) and increased the Weight to
> 0.5 and this folded the peptide into a helix within the 10 ns
> frame I wanted and gave me a PMF. I changed these values
> based on your comment: "...by the time you've added energy
> to one structure metadynamics has moved it back to a similar
> one?" which I didn't quite understand. Isn't it that you add the
> Gaussian to your potential function, wait for newHillFrequency
> steps and add another one? Why should metadynamics move
> it to a similar one? small hills?! Can you please kindly clarify this?
>
> Looking at the end results and the system and applied forces
> of the output, I guess the units are kcal/mol/A (when using
> the RMSD covlar) for the applied forces, right? Do the system
> forces have the same units? Is it normal to observe system forces
> of 5-20 times larger than the biasing (applied?) forces?
> Also, the .pmf file outputs the pmf as a function of the colvar,
> right? I just want to make sure that I got the basics right.
>
> Again, thank you very much for all your time and patience. In
> reporting the other tried and failed results, I thought that (in
> addition to me) maybe someone else may benefit from them
> later.
>
> Cheers,
> Ali
>
>
>
> On Wed, Jun 10, 2009 at 9:57 AM, Giacomo Fiorin <gfiorin_at_seas.upenn.edu>
> wrote:
>>
>> Hi Ali, I see: what you mean (narrowing a1 to an interval of interest) is
>> a
>> good idea, but to do so you also need to set some force constants at
>> the two boundaries. Otherwise, those are just the limits of the grid, and
>> they get expanded over time as the colvar moves outside of them. The
>> fact that lowerBoundary goes negative can be worrisome, but it is only
>> due to the
>> fact that the code managing the grid doesn't know that a1 can't go
>> negative, and allocates a buffer region.
>>
>> Of course, it's never 100% guaranteed that the code is free of bugs,
>> especially at this stage. So, if you do observe negative values of
>> the alpha variable in the .colvars.traj file, I'd like to know.
>>
>> You're applying small hills (0.04 kcal/mol) but very frequently (10
>> steps), so you *are* pumping energy into it. Maybe the helix is just
>> too stiff and by the time you've added energy to one structure
>> metadynamics has moved it back to a similar one?
>>
>> This is why I suggest that you try the harmonic steering: you could,
>> say, set a target of 1 to try and fold the helix completely. Unlike
>> in metadynamics, at least you'd drive it monotonically there, and
>> you'd have to cope with less complications.
>>
>> Otherwise, if your target is the folded helix, and the other
>> structures are not relevant to you, it may simply be happening that
>> because of the system's nature the colvar does not respond as quickly
>> as you want it to, because other degrees of freedom need to be forced
>> (i.e. additional colvars).
>>
>> Ciao
>> Giacomo
>>
>> ---- ----
>> Giacomo Fiorin
>> Center for Molecular Modeling at
>> University of Pennsylvania
>> 231 S 34th Street, Philadelphia, PA 19104-6323
>> phone: (+1)-215-573-4773
>> fax: (+1)-215-573-6233
>> mobile: (+1)-267-324-7676
>> mail: giacomo.fiorin_<at>_gmail.com
>> web: http://www.cmm.upenn.edu/
>> ---- ----
>>
>
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:53:05 CST