[gmx-users] Advice for simulating small DNA

[Joshua Ballanco]

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

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!

- Josh