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
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
