Re: Colvars "tilt" component force constant units for WHAM

From: Robert Elder (rmelder_at_gmail.com)
Date: Fri Jul 22 2011 - 11:06:16 CDT

Hi Giacomo,

Thanks for the response.

On Tue, Jul 12, 2011 at 10:19 PM, Giacomo Fiorin
<giacomo.fiorin_at_gmail.com>wrote:

> Hi Robert, what is numerically speaking, the difference between what you
> get and what you expect?
>

I was seeing that deforming several different molecules was producing
*exactly* the same free energy profile, despite obviously differing profiles
of the colvar I'm using. It ended up being a technical problem with my WHAM
script, and had nothing to do with the force constant, units, or the colvars
module. I just thought I would follow up on this.

> You could change the force constant to another physical unit, for as long
> as the transformation is simply a rescaling, and not a change of variables,
> such as taking the acos(), which has a non-constant derivative. This second
> would require some math, before or after WHAM. So, if you're not OK with
> using the cosine that "tilt" returns, you could try manipulating the code to
> get the tilt angle. Beware though that its gradients are discontinuous
> around 0 and 180, though (as you can visualize easily): the atomic
> gradients of its cosine that you're using, aren't.
>

I decided the easiest thing to do was to use the raw output from colvars
(i.e. the cosine of the angle) and the force constant in the colvars
configuration file for WHAM, and then to transform the x-axis of the free
energy profile from the cosine to the actual angle for easier
interpretation.

>
> Yes, as the manual "suggests", the 1/2 factor is included in the potential.
>

Ah, I was looking at the 2.7b3 manual and couldn't find this. I see where
the 2.8b1 manual very obviously states this (section 10.3.3 for those
following along). Thanks!

> Giacomo
>
>
>
Robert

>
> On Tue, Jul 12, 2011 at 5:48 PM, Robert Elder <rmelder_at_gmail.com> wrote:
>
>> Hi all,
>>
>> I'm a little unclear about the units of the force constant (K) for the
>> new 'tilt' collective variable. The units of the output for this
>> colvar appear to be cos(angle in radians), such that I can extract the
>> angle in degrees using acos(colvars output) * (180/pi).
>>
>> I have been using the raw output from the colvars module (i.e.
>> cos(radians)) and the force constant in my colvars configuration file
>> (no clue what the units are) for WHAM. However, it appears to me that
>> using this force constant may be incorrect, based on a) the expected
>> results of the simulations and b) some testing with lower force
>> constants, which show behavior more closely resembling what I expect.
>> So this made me wonder if I need to convert the K from colvars.conf to
>> some other set of units where the value of K would be lower. Of
>> course, I don't want to choose random numbers until I get desirable
>> results!
>>
>> So, my questions (maybe Jerome or Giacomo can comment):
>>
>> 1. Do I need to convert the tilt force constant to some other set of
>> units for WHAM? Or should I just continue to use the value as it is in
>> my colvars.conf file?
>>
>> 2. Is the factor of 1/2 included in the colvars calculations? That is,
>> is the potential K/2*(x-x0)^2 or simply K*(x-x0)^2? The manual
>> suggests the 1/2 is included, but I'm still unclear.
>>
>> I can provide more details if necessary.
>>
>>
>> Thanks,
>>
>> Robert Elder
>> University of Colorado at Boulder
>> Chemical and Biological Engineering
>>
>>
>

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:20:36 CST