[gmx-users] Advice for simulating small DNA

[Mark Abraham]

Joshua Ballanco wrote:
> 
> On Feb 1, 2009, at 8:22 PM, Mark Abraham wrote:
> 
>> Joshua Ballanco wrote:
>>> On Feb 1, 2009, at 6:48 PM, Mark Abraham wrote:
>>>> Joshua Ballanco wrote:
>>>>> Hi,
>>>>> I'm attempting to model a system involving a small DNA 3-mer. 
>>>>> Without any explicit constraints, the helix begins to come apart 
>>>>> around 0.75 ns to 1 ns into the simulation.
>>>>
>>>> Presumably you have a 3-mer of helices, of which at least one comes 
>>>> apart. Does a single helix in water survive? (Giving a better 
>>>> description of your simulation system would be a good idea!)
>>> Apologies for not being more descriptive... It is a single strand of 
>>> DNA containing 3 A-T base pairs. The system also contains a single 
>>> Arginine residue. Simulating it in water leads to the single DNA 
>>> strand gradually coming apart. With the DNA and water alone, the 
>>> helix stays together much longer, but still eventually comes apart.
>>
>> OK, so your model physics for DNA is intrinsically broken. Where did 
>> you get it?
> 
> The coordinates are from 3DNA. I'm using the terms from the G53a6 force 
> field for DADE and DTHY. As for the H-bond physics, I've thus far been 
> unable to find a good suggestion for how to handle these explicitly 
> using Gromacs. Do you have a recommendation as to where I should be 
> looking? (None of the primary literature I've looked through thus far 
> seems concerned with MD of such short DNA fragments).

No, I've no idea since I don't simulate DNA.

>>>>> 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 restraints.
>>> 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.

>>> 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. :-)

> Thanks again for the pointers. I'm going to try running everything with 
> ffamber to see if it does a better job with the DNA (without the added 
> restraints). I presume the port validated with Gromacs 3.3.1 is still 
> good for 4.0?

Don't know - check out the documentation about the forcefield port - 
search the web.

Mark