[gmx-users] Advice for simulating small DNA
jballanc at gmail.com
Tue Feb 3 02:37:32 CET 2009
On Feb 2, 2009, at 12:26 AM, Mark Abraham wrote:
> No, I've no idea since I don't simulate DNA.
In that case, thank you the help that much more!
>>>>>> So I'm now attempting to add restraints for the base-pair H-
>>>>>> bonds, but I'm having trouble. It seems like no matter what I
>>>>>> try, my system reliably explodes within the first 1 ns. My
>>>>>> constraints look like this:
>>>>>> [ distance_restraints ]
>>>>>> ; ai aj type index type’ low up1 up2 fac
>>>>>> 18 136 1 0 2 0.0 2.0 2.1 1.0
>>>>>> 14 134 1 0 2 0.0 2.0 2.1 1.0
>>>>>> 43 114 1 0 2 0.0 2.0 2.1 1.0
>>>>>> 39 112 1 0 2 0.0 2.0 2.1 1.0
>>>>>> 68 92 1 0 2 0.0 2.0 2.1 1.0
>>>>>> 64 90 1 0 2 0.0 2.0 2.1 1.0
>>>>>> I've tried pre-equilibrating for up to 100 ps, but even that
>>>>>> doesn't prevent the system from eventually exploding.
>>>>> Your .mdp settings for distance restraints may also be relevant
>>>>> here - not least in setting the existence and magnitude of these
>>>> As I understand, the only relevant lines are:
>>>> constraints = all-bonds
>>>> integrator = md
>>>> disre = simple
>>> disre-fc and others are also relevant. See manual chapter 7.
>> Thanks for the pointer. I had overlooked most of the options there,
>> since I'm not actually doing anything related to NMR. (That'll
>> teach me to read more carefully!) Unfortunately, playing around
>> with this, disre-tau, disre-weighting, and the weighting factors
>> for each bond have not, so far, avoided the explosion.
> OK, that's no longer surprising - distance restraints will not
> usefully fix a broken model physics.
Well, yes, but I also wouldn't expect them to break the broken physics
further... I realize the system I was using originally was rather
unphysical, but the DNA helix at least was at least *mostly* holding
together. When I add the distance restraints, even with very large
multipliers, the seem to serve only to tear apart the helix. Odd...
>>>> For PME I was using:
>>>> coulombtype = PME
>>>> rlist = 0.55
>>>> rcoulomb = 0.55
>>>> rvdw = 0.55
>>>> fourierspacing = 0.1375
>>> I agree with Justin that these are very weird for normal usage.
>> Thanks for pointing that out. I'm relatively new with Gromacs, and
>> hastily reduced these values to fix the relatively small box my
>> system fits in. I doubled the short box dimension (triclinic; was --
>> > 2.0, 2.1, 1.1 now --> 2.0, 2.1, 2.0) and increased the radii to
>> the (as far as I can tell) more recommended values:
>> coulombtype = PME
>> rlist = 0.9
>> rcoulomb = 0.9
>> rvdw = 0.9
>> fourierspacing = 0.12
> Well, that's more like it. Values for these parameters are actually
> intrinsic to the forcefield parametrization process, and one should
> vary them only with caution. This algorithmic constraint sets a
> minimum size for the simulation, of course.
> When using PBC, just fitting your system into a box doesn't address
> the real issue. In a real solution this 3-mer would be close to
> infinite dilution, which can't be modeled without a serious chunk of
> solvent around it. This consideration dwarfs the algorithmic one I
> refer to above.
>> Unfortunately, even with all of these changes, I'm still getting an
>> explosion (and my simulation is quite a bit slower).
> Anyone can get quick random numbers - you don't even need a
> simulation package :-P There's no substitute for background
> literature reading, doing tutorials, and experimenting with
> preparedness for failure. :-)
Agreed, and thank you for all the help. I've decided to get a bit
creative, and work around the whole short-DNA fragment issue all
together. The new approach will also reduce the number of simulations
I'm looking at doing from ~6000 to ~20, so I don't have to worry about
the larger boxes taking so much longer to simulate.
Again, thanks for everything!
More information about the gromacs.org_gmx-users