From: Michelle Kuttel (mkuttel_at_cs.uct.ac.za)
Date: Thu Nov 24 2011 - 07:24:58 CST
Hi Jerome
I do know that relaxed maps are complicated to calculate correctly (one requirement being that all other degrees of freedom have to be in their most favourable conformation). However, correctness aside, I really still don't understand the problem behind my original question, which was "Why does changing the TargetNumStages, or TargetNumSteps, require a higher force to produce the same profile?" (whether the profile is 'correct' or not is a separate concern). This behaviour does not make sense to me and it really is a question about how the method is programmed in NAMD. However, this discussion has reached the stage where it is now probably only of interest to me .... ;)
Many thanks for your time.
Michelle
-----------------------------
Michelle Kuttel
mkuttel_at_cs.uct.ac.za
PH: +27 21 6505107
-----------------------------
On 24 Nov 2011, at 2:33 PM, Jérôme Hénin wrote:
> Dear Michelle,
>
> After inspecting your files, I think you are seeing an intrinsic
> limitation of relaxed scans, that is, energy minimization algorithms
> get trapped easily. In this case, the intramolecular H-bonds form or
> not, depending on the precise orientation of the approaching OH group.
> Getting the relaxed potential energy surface requires a lot of
> attention and possibly manual adjustments to remain in the
> (constrained) global minimum. As you have seen, MD approaches at room
> temperature do not suffer from this quite as much.
>
> A few points about your setup: instead of defining the dihedral based
> on a hydrogen atom, I would tend to choose the hydroxyl oxygen, which
> is more directly representative of the conformational change.
>
> More importantly, the value you used as a starting point (-180) is far
> from the starting configuration, so that the first stage of
> minimization required extensive rearrangement. Instead, I would scan
> from -60 to 300 (using the H atom definition, or -300 to 60 if using
> the oxygen). As a check, you can also scan from -300 to 420, and see
> if both full turns produce the same results. I also tried backwards,
> which produces different results (H-bonds have a greater tendency to
> break).
>
> The bottom line is, if ABF produce the result you need, then use it by
> all means.
>
> Best,
> Jerome
>
> On 24 November 2011 10:51, Jérôme Hénin <jhenin_at_ifr88.cnrs-mrs.fr> wrote:
>> Dear Michelle:
>>
>> Maybe I should have a closer look at this. Can you send me off-list
>> the files I need to reproduce your torsional scan?
>>
>> Thanks,
>> Jerome
>>
>> On 24 November 2011 06:33, Michelle Kuttel <mkuttel_at_cs.uct.ac.za> wrote:
>>> Dear Aron
>>>
>>> Hi Michelle,
>>>
>>> From what you are saying am I to understand that your glucose is just in a
>>> vacuum? So there shouldn't be any bad interactions with any other
>>> molecules?
>>>
>>> Yes, that is right!
>>>
>>> If so, the only reason I could think of would be that your colvar is not
>>> defined properly (e.g. the groups are in the wrong order, or the atom
>>> numbers don't match), but you say you've done ABF with exactly the same
>>> colvar definition and it worked fine?
>>>
>>> Yes, worked beautifully.
>>>
>>> If so then I'd be inclined to agree with you. Sadly I've never used the
>>> minimization routine to calculate anything, so I can't be more helpful.
>>>
>>> Well, it is actually helpful that you agree (i.e. I'm not being a complete
>>> idiot). ;)
>>> Thanks
>>> Michelle
>>>
>>> ~Aron
>>>
>>> On Wed, Nov 23, 2011 at 3:03 PM, Michelle Kuttel <mkuttel_at_cs.uct.ac.za>
>>> wrote:
>>>>
>>>> Hi all
>>>>
>>>>
>>>> On 23 Nov 2011, at 6:41 PM, Jérôme Hénin wrote:
>>>>
>>>>> Hi guys,
>>>>>
>>>>> On 23 November 2011 16:28, Giacomo Fiorin <giacomo.fiorin_at_gmail.com>
>>>>> wrote:
>>>>>> Hi Michelle. What you describe seems to me that the angle gets stuck
>>>>>> into a
>>>>>> local minimum. This is pretty typical when you minimize. Then when
>>>>>> you try
>>>>>> the next window, you don't get the results you want because your angle
>>>>>> is
>>>>>> stuck.
>>>>>>
>>>>>> I would start all windows from the same initial structure and change
>>>>>> the
>>>>>> center of the static restraint, rather than changing it on the fly
>>>>>> while you
>>>>>> minimize. That would be most appropriate for dynamics.
>>>>>
>>>>> I am not sure that will work so well. The problem is more fundamental:
>>>>> if the potential is corrugated enough that you get stuck in local
>>>>> minima, then the relaxed scan is always tricky. I run into this once,
>>>>> and one test I did was to run the scan both ways to see if I sampled
>>>>> the same conformations - and if not, which one had lower energy.
>>>>>
>>>>> You could try scanning this angle while doing MD, to give the molecule
>>>>> a chance to optimize itself beyond the nearest local minima.
>>>>>
>>>>> Once the "true" minima are identified, you might have to run the scan
>>>>> in pieces, starting from the properly minimized structures and
>>>>> sampling the barriers around going both ways. In the worst case, the
>>>>> landscape is too complex to be sampled like this, then you'd need some
>>>>> kind of MC search rather than just a torsional scan.
>>>>
>>>> Actually, this is a very simple dihedral angle in a very small molecule:
>>>> the primary alcohol in glucose. I was running it as a test, to validate
>>>> the force field in NAMD (it's a CHARMM ff) and procedure before I continue
>>>> with more complex calculations. It really is one of the simplest cases you
>>>> can get, so colvars with "harmonic" *should* work for calculating a relaxed
>>>> map. BTW, the ABF routine works beautifully on the same dihedral, so I
>>>> think that I am forced to conclude that "harmonic" does not actually work
>>>> for calculating relaxation maps. Or at least not reliably.
>>>>
>>>>>
>>>>>
>>>>>
>>>>>> On Wed, Nov 23, 2011 at 9:16 AM, Michelle Kuttel <mkuttel_at_cs.uct.ac.za>
>>>>>> wrote:
>>>>>>>
>>>>>>> The plot thickens....
>>>>>>> After much experimenting, I have finally determined that, for this
>>>>>>> particular dihedral, the optimal "harmonic" configuration is:
>>>>>>> harmonic {
>>>>>>> name scanOmega
>>>>>>> colvars Omega
>>>>>>> centers -180
>>>>>>> targetCenters 180
>>>>>>> targetNumSteps 500 #these
>>>>>>> targetNumStages 72 #are the values
>>>>>>> forceConstant 0.005 # in question
>>>>>>> }
>>>>>>> So ,it seems that my previous values were far too "high". However, I
>>>>>>> still don't understand why, if I increase targetNumStages and
>>>>>>> targetNumSteps, I need to increase the force constant as well.
>>>>>>> "targetNumSteps" should just be number of minimization steps at a
>>>>>>> particular force, surely, and not affect the applied force, right?
>>>>>>> But it
>>>>>>> does!
>>>>>
>>>>> Note that the force constant is in unit of kcal/mol/square colvar
>>>>> unit, where the effective colvar unit is given by the width parameter.
>>>>> Here, it is 0.5 degrees, which is very small, meaning that a force
>>>>> constant of 100 is really huge (deviating from the restraint center by
>>>>> 0.5 degrees costs 100 kcal/mol of energy). A looser restraint leaves
>>>>> more leeway for various degrees of freedom to relax, possibly yielding
>>>>> a better optimization.
>>>>>
>>>>> The reason targetNumStages matters is, you may not get trapped in the
>>>>> same local minima depending on the size of the "hops" in the reference
>>>>> torsional angle.
>>>>
>>>> But this is a simple dihedral, with only three minima (at -60, 60 and 180
>>>> degrees) - and these minima are not where it gets stuck....
>>>> I am still not clear on why "targetNumSteps" affects the force....?
>>>>
>>>> Many thanks
>>>> Michelle
>>>>
>>>>>
>>>>> Best,
>>>>> Jerome
>>>>>
>>>>>>> The use of "harmonic" to create relaxed energy maps is not documented
>>>>>>> in
>>>>>>> the NAMD User Guide, other than the comment:
>>>>>>> "The harmonic biasing method may be used to enforce fixed or moving
>>>>>>> restraints, including variants of Steered and Targeted MD. Within
>>>>>>> energy
>>>>>>> minimization runs, it allows for restrained minimization, e.g. to
>>>>>>> calculate
>>>>>>> relaxed potential energy surfaces.".
>>>>>>> So, I would really appreciate it if someone can explain this to me.
>>>>>>> Pretty please?
>>>>>>> Michelle
>>>>>>>
>>>>>>> On 21 Nov 2011, at 12:45 PM, Michelle Kuttel wrote:
>>>>>>>
>>>>>>> Hello
>>>>>>> I am trying to calculate a simple relaxed energy profile for a
>>>>>>> dihedral
>>>>>>> angle (this is a simple sanity check - I can compare the profile to
>>>>>>> one I
>>>>>>> calculated some time ago in CHARMM).
>>>>>>> I am using colvars and harmonic restraints to do this. My
>>>>>>> colvarsConfig
>>>>>>> file looks like this:
>>>>>>> ------------------------------------
>>>>>>> colvarsTrajFrequency 10000
>>>>>>> colvar {
>>>>>>> name Omega
>>>>>>> width 0.5
>>>>>>> dihedral {
>>>>>>> group1 {
>>>>>>> atomnumbers { 7 }
>>>>>>> }
>>>>>>> group2 {
>>>>>>> atomnumbers { 5 }
>>>>>>> }
>>>>>>> group3 {
>>>>>>> atomnumbers { 20 }
>>>>>>> }
>>>>>>> group4 {
>>>>>>> atomnumbers { 21 }
>>>>>>> }
>>>>>>> }
>>>>>>>
>>>>>>> }
>>>>>>> harmonic {
>>>>>>> name scanOmega
>>>>>>> colvars Omega
>>>>>>> centers -180
>>>>>>> targetCenters 180
>>>>>>> targetNumSteps 10000
>>>>>>> targetNumStages 36
>>>>>>> forceConstant 100
>>>>>>> }
>>>>>>> ----------------------------
>>>>>>> In my main file, I then have the line:
>>>>>>> minimize 370000
>>>>>>> This works OK (though the energies for the first few stages are a bit
>>>>>>> high). HOWEVER, when I change "targetNumStages 36" to
>>>>>>> "targetNumStages 72"
>>>>>>> and "minimize 370000" to "minimize 730000", suddenly I need a much
>>>>>>> higher
>>>>>>> force of 200 for it to work. If I don't increase the force, the omega
>>>>>>> angle
>>>>>>> (recorded in the .traj file) "stalls" at -145 degrees. The angle
>>>>>>> does not
>>>>>>> rotate further than this on subsequent stages.
>>>>>>> Why would increasing the number of stages require a higher force? I
>>>>>>> thought that the width of the omega colvar would still be the same
>>>>>>> (0.5
>>>>>>> degrees) and therefore the scaling should be the same?
>>>>>>> I think that I do not understand the way forces are applied under
>>>>>>> harmonic
>>>>>>> restraints and I would really appreciate some insight!
>>>>>>> Many thanks
>>>>>>> Michelle
>>>>>>> ---------------------------------------------
>>>>>>> Dr Michelle Kuttel
>>>>>>> Department of Computer Science
>>>>>>> University of Cape Town
>>>>>>> Cape Town
>>>>>>> South Africa
>>>>>>> mkuttel_at_cs.uct.ac.za
>>>>>>> PH: +27 21 6505107
>>>>>>> ---------------------------------------------
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>
>
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:21:02 CST