Re: Scaled forces, colvars and harmonic con/restraints

From: Jérôme Hénin (jhenin_at_ifr88.cnrs-mrs.fr)
Date: Thu Nov 24 2011 - 03:51:51 CST

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:01 CST