[gmx-users] Advice for simulating small DNA

[Justin A. Lemkul]

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

Read the primary literature references for Gromos96 53a6; I believe the rvdw 
parameter should be 1.4 nm to keep consistent with the parameterization scheme.

Are your box dimensions really adequate?  A standard DNA helix should be about 2 
nm wide, so your system may be seeing its periodic images if your box is only 
about 2.0 nm across its shortest dimension.  Use, for example, editconf -c -d 
1.0 to generate sufficient box dimensions.

-Justin

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

-- 
========================================

Justin A. Lemkul
Graduate Research Assistant
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================