# Re: Scaled forces, colvars and harmonic con/restraints

From: Michelle Kuttel (mkuttel_at_cs.uct.ac.za)
Date: Thu Nov 24 2011 - 09:39:56 CST

Dear Jerome

Thank you! This now makes sense, which is a relief. The loop you suggested works, except for the slight difference that the collective variable values are now printed twice to the .colvars.traj. file.
Many thanks for paying this rather obscure question such attention.
It is much appreciated.
regards
Michelle

-----------------------------
Michelle Kuttel
mkuttel_at_cs.uct.ac.za
PH: +27 21 6505107
-----------------------------

On 24 Nov 2011, at 4:18 PM, Jérôme Hénin wrote:

> Hi Michelle:
>
> I have now (finally) reproduced the behavior that prompted your
> question in the first place. What I see is that things work as long as
> the energy is not too well-minimized. Past a certain number of
> minimization steps, when the target dihedral changes, the minimizer
> gets stuck and oscillates. The reason for this, I've finally figured
> out, is that the minimizer doesn't adapt well to sudden changes in the
> energy function.
>
> Fortunately, there is a simple solution to this (which I should have
> thought of in the first place). Reset the minimizer at each stage by
> calling it like this:
>
> set numStages 36
> set numSteps 1000
>
> for {set i 0} {\$i <= \$numStages} {incr i} {
> minimize \$numSteps
> }
>
> This should work every time (notwithstanding the problem of local
> minima). Please let me know if it doesn't.
>
> Thank you for reporting this problem.
> Jerome
>
>
>
> On 24 November 2011 14:24, Michelle Kuttel <mkuttel_at_cs.uct.ac.za> wrote:
>> 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.
>>>
>>> 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.
>>>>>>>>> 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 : Wed Feb 29 2012 - 15:57:58 CST