[gmx-users] Error in DNA Simulation

[Justin Lemkul]

On 11/4/12 9:15 AM, Mehdi Bagherpour wrote:
> Dear all,
>
> I used Gromacs to Simulate DNA(15bp)  with CHARMM27 force-field in 50ns .My
> sistem has about 10000 water molecules and  outside of the DNA and the edge
> of the box is 1nm. At the end of simulation the DNA is destroyed.(e.g. DNA
> has bended  and last base pair has separated from helix)
> I used 8 processor system for simulation.

You have a large number of incorrect settings in your .mdp file.  I will only 
comment on the md.mdp file, but the same apply to NVT and NPT runs.

<snip>

> Ubq  md.mdp
>
> ; 7.3.3 Run Control
> integrator              = md                    ; md integrator
> tinit                   = 0                     ; [ps] starting time for run
> dt                      = 0.002                 ; [ps] time step for
> integration
> nsteps                  = 2500000               ; maximum number of steps
> to integrate, 0.002 * 14,000,000 = 28,000 ps
> comm_mode               = Linear                ; remove center of mass
> translation
> nstcomm                 = 1                     ; [steps] frequency of mass
> motion removal
> comm_grps               = DNA     Sol       Ion        MG   ; group(s) for
> center of mass motion removal
>

You should not remove COM motion of individual components.  Set this to "System."

> ; 7.3.8 Output Control
> nstxout                 = 500       ; [steps] freq to write coordinates to
> trajectory
> nstvout                 = 500       ; [steps] freq to write velocities to
> trajectory
> nstfout                 = 500       ; [steps] freq to write forces to
> trajectory
> nstlog                  = 100           ; [steps] freq to write energies to
> log file
> nstenergy               = 500           ; [steps] freq to write energies to
> energy file
> nstxtcout               = 500           ; [steps] freq to write coordinates
> to xtc trajectory
> xtc_precision           = 1000          ; [real] precision to write xtc
> trajectory
> xtc_grps                = System        ; group(s) to write to xtc
> trajectory
> energygrps              = System        ; group(s) to write to energy file
>
> ; 7.3.9 Neighbor Searching
> nstlist                 = 1             ; [steps] freq to update neighbor
> list

You lose a lot of performance doing this.  You only need nstlist = 1 during EM. 
  It won't cause a crash, but your simulations will be unnecessarily slow.  A 
value of 5 or 10 is sufficient.

> ns_type                 = grid          ; method of updating neighbor list
> pbc                     = xyz           ; periodic boundary conditions in
> all directions
> rlist                   = 1           ; [nm] cut-off distance for the
> short-range neighbor list
>
> ; 7.3.10 Electrostatics
> coulombtype             = PME           ; Particle-Mesh Ewald electrostatics
> rcoulomb                = 1           ; [nm] distance for Coulomb cut-off
>
> ; 7.3.11 VdW
> vdwtype                 = cut-off       ; twin-range cut-off with rlist
> where rvdw >= rlist
> rvdw                    = 1           ; [nm] distance for LJ cut-off
> DispCorr                = EnerPres      ; apply long range dispersion
> corrections for energy
>

Your cutoff settings are all wrong for CHARMM27.  You need:

rlist = 1.2
rlistlong = 1.4
vdw_type = switch
rvdw_switch = 0.8
rvdw = 1.2
rcoulomb = 1.2

> ; 7.3.13 Ewald
> fourierspacing          = 0.12          ; [nm] grid spacing for FFT grid
> when using PME
> pme_order               = 4             ; interpolation order for PME, 4 =
> cubic
> ewald_rtol              = 1e-5          ; relative strength of
> Ewald-shifted potential at rcoulomb
>
> ; 7.3.14 Temperature Coupling
> tcoupl                  = v-rescale                   ; temperature
> coupling with Nose-Hoover ensemble
> tc_grps                 = DNA        Sol       Ion
> MG                 ; groups to couple seperately to temperature bath
> tau_t                   = 0.1       0.1       0.1
> 0.1                   ; [ps] time constant for coupling
> ref_t                   = 300       300       300
> 300                  ; [K] reference temperature for coupling
>

You should not couple ions separately from the solvent.

http://www.gromacs.org/Documentation/Terminology/Thermostats

> ; 7.3.15 Pressure Coupling
> pcoupl                  = parrinello-rahman     ; pressure coupling where
> box vectors are variable
> pcoupltype              = isotropic             ; pressure coupling in
> x-y-z directions
> tau_p                   = 2.0                   ; [ps] time constant for
> coupling
> compressibility         = 4.5e-5                ; [bar^-1] compressibility
> ref_p                   = 1.0                   ; [bar] reference pressure
> for coupling
>
> ; 7.3.17 Velocity Generation
> gen_vel                 = no            ; velocity generation turned off
>
> ; 7.3.18 Bonds
> constraints             = hbonds     ; convert all bonds to constraints
> constraint_algorithm    = LINCS         ; LINear Constraint Solver
> continuation            = yes           ; apply constraints to the start
> configuration
> lincs_order             = 4             ; highest order in the expansion of
> the contraint coupling matrix
> lincs_iter              = 1             ; number of iterations to correct
> for rotational lengthening
> lincs_warnangle         = 30            ; [degrees] maximum angle that a
> bond can rotate before LINCS will complain
> *************************************************************************
> Also when my run is crashed,i used from followed command:
>
> tpbconv -s wi_md.tpr -extend 40000 -o wi_md_1.tpr
> mdrun -s wi_md_1.tpr -cpi wi_md.cpt -deffnm wi_md_1 -append yes
> *************************************************************************
>

If the run crashes (presumably with LINCS warnings or something else that 
indicates the system is blowing up), it is usually futile to try to restart it. 
  There is something unsound about the system.  Given the extent of the above 
errors, it seems very likely that you're simply breaking the force field.

-Justin

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

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

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