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

[ms]

Hi,

It's a bit long but bear with me if you can. I'm getting quite mad with 
this. Thanks :)

Now. I am using Gromacs (4.0.7 currently) to create a custom 
coarse-grained, no-solvent MD peptide model -one that contains the 
backbone but has only C-alphas, no side chains.

To enforce chirality in such a toy, I thought that a simple (naive?) but 
functional idea could be that of enforcing a not-too-hard and 
wide-bottomed dihedral restrain on the phi angle, like that:

[ dihedral_restraints ]
; ai   aj    ak    al  type  label  phi  dphi  kfac  power
; phi C'(n-1) - N - CA - C'
   3     5       7       8       1        1      -90     75    2      1

so that the system is free to move in a large phi-angle space without 
problems, but cannot go in the forbidden, symmetrical glycine-like 
space. In testing, I set the force constant small (10) so that the 
chain, even if with some difficulty, should be able to bend and move by 
spending some energy in the "forbidden" zones, if needed to achieve 
structure, etc.

Now, the point is that *when* it works, it works quite well. But 
sometimes -with "sometimes" being from 0% to 90% of times, depending on 
details of what I am simulating- it totally doesn't. That is: sometimes, 
starting from an extended structure, the chain simply freezes. It is not 
static: it wiggles a bit here and there, but it stays mostly extended. 
Energies plots (g_energy output xvg's) are on average flat: that is, 
they show the right noise for all components, but no structure: nothing 
happens, for as much as 50 ns.

Sometimes instead it works, and then you see that everything works 
alright; the peptide structures, potential energy goes down and down, 
etc.etc.

There are no differences in minimization energies etc. from a simulation 
that freezes and one that works. Moving the restraint center doesn't 
help much. I *know* it's the restraints because if I keep everything the 
same but I lower the force constant or remove restraints altogether, the 
freezing bug happens no more -but then I have little or no more 
chirality, of course. I have the following further facts:
- It tends to happen more often for longer chains: 10 or 20-residue long 
chains work often OK, but a 40-long chain is difficult to get "right".
- It *never* happens, to my knowledge, if I put more than 1 chain to 
cluster in the box.
- The Ramachandran plot of the extended frozen thing is "right", meaning 
that it populates a reasonable (mostly beta) angle space. During 
trajectory angles move quite a bit but stay in the right (top-left) zone 
of the plot.
- "Good"-behaving simulation often start with a short "frozen" stretch: 
then the dihedral restraint energy suddenly jumps to some above noise 
positive value, to accomodate some bending I guess, and all goes 
downhill from there.
- In frozen simulations, there's almost no LJ energy: the chain does not 
self-interact.

All of this, added to the apparent "randomness" of the pathology, makes 
me think that the chain simply is stuck in a situation where it has a 
too high barrier to properly bend -a barrier which is not a "soft" one 
as the harmonic potential of the restrain would make me think. This 
would explain why putting more peptides and letting them bind together 
avoids the problem: the potential energy of binding compensates for the 
bending when you have more than one thing; but sometimes a single chain 
cannot find itself to interact with, because it should bend correctly to 
do so.

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:
- How can I try to understand in full what is going on
- If my diagnosis sounds right, how can I start to correct the problem.
- If the problem is my naive approach to chirality, any suggestion is 
welcome.

Thanks A LOT if you've been patient enough to read until here! And 
thanks even *more* if you can help!

Yours truly,
Massimo