[gmx-users] The trouble with dihedral restraints: frozen peptide backbones

[ms]

Ok, I can update this situation (see below quote):

On 18/10/10 03:30, Mark Abraham wrote:
>> And then, I read on the Gromacs website *this*
>> http://www.gromacs.org/Documentation/How-tos/Dihedral_Restraints
>> " 2. The manual is a bit unclear about whether this type of
>> dihedral restraint is stable for use near 180 degrees. Chris
>> Neale has found that everything appears to behave normally and
>> as expected over the entire range of dihedral angles including
>> 180 degrees. However, one must avoid the situation in which the
>> actual dihedral is close to 180deg away from the restrained dihedral."
>>
>> And I know that somehow the potential is discontinuous and/or
>> anyway not intended to be well-behaved everywhere:
>> http://www.mail-archive.com/gmx-users@gromacs.org/msg12353.html
>>
>> ">  On Tue, Feb 26, 2008 at 11:37 AM, David Mobley<[EMAIL
>> PROTECTED]>  wrote:
>>>>   The other way of putting what Mark said is that phi is
>> only meaningful on some range (-pi to pi, or 0 to 2pi, depending
>> on how you define it)
>>>>   and so what you require is that the potential be
>> harmonic for the
>>>>   region in which phi is meaningful. You don't care what
>> happens outside that. Or, in this case, you handle the issue by
>> mapping phi values  outside the allowed range back into
>> that allowed range."
>>
>> (I don't get what does it mean to "map back" phi values, by the way).
>>
>> I am no MD expert, so maybe my hypothesis is totally off, but I
>> am reasonably sure that either I do not understand correctly how
>> dihedral restraints work, or I am suffering by the discontinuity
>> of the way the potential is implemented (Or both). I would like
>> to know:

I managed to patch the dihedral restraining part so to mitigate the 
issue. The problem lied in the 180-degrees discontinuity in the force, 
where it jumps from a highly positive to a highly negative value, like 
(Y=force, X=angle):

    /|
   / |
  /  |  /
/   | /
     |/

I have cut the ends close to the phi=180 and "connected" the force 
values like that:

   /\
  /  \    /
/    \  /
       \/

and then redefined the energy to be the integral of that line.

This way it works -you can have decent restraining along all phi angles 
without risking to freeze, provided the energy barrier at 180 is not too 
high.

The patch looks really ugly -it has been kludged in quite a hurry and... 
well, looks obscene, but if someone wants to look at it, let me know!

m.

-- 
Massimo Sandal, Ph.D.
http://devicerandom.org