[gmx-users] Advice for simulating small DNA

[Mark Abraham]

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?

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

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

Mark