<br><br>----- Original Message -----<br>From: ms &lt;devicerandom@gmail.com&gt;<br>Date: Monday, October 18, 2010 8:23<br>Subject: [gmx-users] The trouble with dihedral restraints: frozen peptide backbones<br>To: Discussion list for GROMACS users &lt;gmx-users@gromacs.org&gt;<br><br>&gt; Hi,<br>&gt; <br>&gt; It's a bit long but bear with me if you can. I'm getting quite <br>&gt; mad with this. Thanks :)<br>&gt; <br>&gt; Now. I am using Gromacs (4.0.7 currently) to create a custom <br>&gt; coarse-grained, no-solvent MD peptide model -one that contains <br>&gt; the backbone but has only C-alphas, no side chains.<br>&gt; <br>&gt; To enforce chirality in such a toy, I thought that a simple <br>&gt; (naive?) but functional idea could be that of enforcing a not-<br>&gt; too-hard and wide-bottomed dihedral restrain on the phi angle, <br>&gt; like that:<br>&gt; <br>&gt; [ dihedral_restraints ]<br>&gt; ; ai&nbsp;&nbsp; aj&nbsp;&nbsp;&nbsp; ak&nbsp;&nbsp;&nbsp; <br>&gt; al&nbsp; type&nbsp; label&nbsp; phi&nbsp; dphi&nbsp; kfac&nbsp; power<br>&gt; ; phi C'(n-1) - N - CA - C'<br>&gt; &nbsp; 3&nbsp;&nbsp;&nbsp;&nbsp; <br>&gt; 5&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 7&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 8&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; -90&nbsp;&nbsp;&nbsp;&nbsp; 75&nbsp;&nbsp;&nbsp; 2&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 1<br><br>Try mdrun -debug on a short run. What does it have to say about the parameter values for this dihedral restraint? Likewise gmxdump -s on the .tpr (though you will have to search around for F_DIHRES among other things, probably).<br><br>Mark<br><br>&gt; <br>&gt; so that the system is free to move in a large phi-angle space <br>&gt; without problems, but cannot go in the forbidden, symmetrical <br>&gt; glycine-like space. In testing, I set the force constant small <br>&gt; (10) so that the chain, even if with some difficulty, should be <br>&gt; able to bend and move by spending some energy in the "forbidden" <br>&gt; zones, if needed to achieve structure, etc.<br>&gt; <br>&gt; Now, the point is that *when* it works, it works quite well. But <br>&gt; sometimes -with "sometimes" being from 0% to 90% of times, <br>&gt; depending on details of what I am simulating- it totally <br>&gt; doesn't. That is: sometimes, starting from an extended <br>&gt; structure, the chain simply freezes. It is not static: it <br>&gt; wiggles a bit here and there, but it stays mostly extended. <br>&gt; Energies plots (g_energy output xvg's) are on average flat: that <br>&gt; is, they show the right noise for all components, but no <br>&gt; structure: nothing happens, for as much as 50 ns.<br>&gt; <br>&gt; Sometimes instead it works, and then you see that everything <br>&gt; works alright; the peptide structures, potential energy goes <br>&gt; down and down, etc.etc.<br>&gt; <br>&gt; There are no differences in minimization energies etc. from a <br>&gt; simulation that freezes and one that works. Moving the restraint <br>&gt; center doesn't help much. I *know* it's the restraints because <br>&gt; if I keep everything the same but I lower the force constant or <br>&gt; remove restraints altogether, the freezing bug happens no more -<br>&gt; but then I have little or no more chirality, of course. I have <br>&gt; the following further facts:<br>&gt; - It tends to happen more often for longer chains: 10 or 20-<br>&gt; residue long chains work often OK, but a 40-long chain is <br>&gt; difficult to get "right".<br>&gt; - It *never* happens, to my knowledge, if I put more than 1 <br>&gt; chain to cluster in the box.<br>&gt; - The Ramachandran plot of the extended frozen thing is "right", <br>&gt; meaning that it populates a reasonable (mostly beta) angle <br>&gt; space. During trajectory angles move quite a bit but stay in the <br>&gt; right (top-left) zone of the plot.<br>&gt; - "Good"-behaving simulation often start with a short "frozen" <br>&gt; stretch: then the dihedral restraint energy suddenly jumps to <br>&gt; some above noise positive value, to accomodate some bending I <br>&gt; guess, and all goes downhill from there.<br>&gt; - In frozen simulations, there's almost no LJ energy: the chain <br>&gt; does not self-interact.<br>&gt; <br>&gt; All of this, added to the apparent "randomness" of the <br>&gt; pathology, makes me think that the chain simply is stuck in a <br>&gt; situation where it has a too high barrier to properly bend -a <br>&gt; barrier which is not a "soft" one as the harmonic potential of <br>&gt; the restrain would make me think. This would explain why putting <br>&gt; more peptides and letting them bind together avoids the problem: <br>&gt; the potential energy of binding compensates for the bending when <br>&gt; you have more than one thing; but sometimes a single chain <br>&gt; cannot find itself to interact with, because it should bend <br>&gt; correctly to do so.<br>&gt; <br>&gt; And then, I read on the Gromacs website *this*<br>&gt; http://www.gromacs.org/Documentation/How-tos/Dihedral_Restraints<br>&gt; " 2. The manual is a bit unclear about whether this type of <br>&gt; dihedral restraint is stable for use near 180 degrees. Chris <br>&gt; Neale has found that everything appears to behave normally and <br>&gt; as expected over the entire range of dihedral angles including <br>&gt; 180 degrees. However, one must avoid the situation in which the <br>&gt; actual dihedral is close to 180deg away from the restrained dihedral."<br>&gt; <br>&gt; And I know that somehow the potential is discontinuous and/or <br>&gt; anyway not intended to be well-behaved everywhere:<br>&gt; http://www.mail-archive.com/gmx-users@gromacs.org/msg12353.html<br>&gt; <br>&gt; "&gt; On Tue, Feb 26, 2008 at 11:37 AM, David Mobley &lt;[EMAIL <br>&gt; PROTECTED]&gt; wrote:<br>&gt; &gt; &gt;&nbsp; The other way of putting what Mark said is that phi is <br>&gt; only meaningful on some range (-pi to pi, or 0 to 2pi, depending <br>&gt; on how you define it)<br>&gt; &gt; &gt;&nbsp; and so what you require is that the potential be <br>&gt; harmonic for the<br>&gt; &gt; &gt;&nbsp; region in which phi is meaningful. You don't care what <br>&gt; happens outside that. Or, in this case, you handle the issue by <br>&gt; mapping phi values&nbsp; outside the allowed range back into <br>&gt; that allowed range."<br>&gt; <br>&gt; (I don't get what does it mean to "map back" phi values, by the way).<br>&gt; <br>&gt; I am no MD expert, so maybe my hypothesis is totally off, but I <br>&gt; am reasonably sure that either I do not understand correctly how <br>&gt; dihedral restraints work, or I am suffering by the discontinuity <br>&gt; of the way the potential is implemented (Or both). I would like <br>&gt; to know:<br>&gt; - How can I try to understand in full what is going on<br>&gt; - If my diagnosis sounds right, how can I start to correct the <br>&gt; problem.- If the problem is my naive approach to chirality, any <br>&gt; suggestion is welcome.<br>&gt; <br>&gt; Thanks A LOT if you've been patient enough to read until here! <br>&gt; And thanks even *more* if you can help!<br>&gt; <br>&gt; Yours truly,<br>&gt; Massimo<br>&gt; -- <br>&gt; gmx-users mailing list&nbsp;&nbsp;&nbsp; gmx-users@gromacs.org<br>&gt; http://lists.gromacs.org/mailman/listinfo/gmx-users<br>&gt; Please search the archive at <br>&gt; http://www.gromacs.org/Support/Mailing_Lists/Search before posting!<br>&gt; Please don't post (un)subscribe requests to the list. Use the <br>&gt; www interface or send it to gmx-users-request@gromacs.org.<br>&gt; Can't post? Read http://www.gromacs.org/Support/Mailing_Lists