[gmx-users] Advice for simulating small DNA

[Joshua Ballanco]

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

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

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

Unfortunately, even with all of these changes, I'm still getting an  
explosion (and my simulation is quite a bit slower).

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?

Thanks!

- Josh