[gmx-users] Error in DNA Simulation

Mehdi Bagherpour mbagherpour7749 at gmail.com
Sun Nov 4 15:15:06 CET 2012


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.
My mdp files for NVT,NPT and MD files are:
*********************************************************************************
Ubq nvt.mdp

; 7.3.2 Preprocessing
define                  = -DPOSRES      ; defines to pass to the
preprocessor

; 7.3.3 Run Controlnstcgsteep
integrator              = md                    ; md integrator
tinit                   = 0                     ; [ps] starting time for run
dt                      = 0.002                 ; [ps] time step for
integration
nsteps                  = 50000                 ; maximum number of steps
to integrate, 0.002 * 25,000 = 50 ps
comm_mode               = Linear                ; remove center of mass
translation
nstcomm                 = 1                     ; [steps] frequency of mass
motion removal

; 7.3.8 Output Control
nstxout                 = 25000         ; [steps] freq to write coordinates
to trajectory
nstvout                 = 25000         ; [steps] freq to write velocities
to trajectory
nstfout                 = 25000         ; [steps] freq to write forces to
trajectory
nstlog                  = 100           ; [steps] freq to write energies to
log file
nstenergy               = 100           ; [steps] freq to write energies to
energy file
nstxtcout               = 100           ; [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
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

; 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 Berendsen-thermostat
tc_grps                 = DNA      Ion      Sol      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

; 7.3.17 Velocity Generation
gen_vel                 = yes           ; generate velocities according to
Maxwell distribution of temperature
gen_temp                = 300           ; [K] temperature for Maxwell
distribution
gen_seed                = -1            ; [integer] used to initialize
random generator for random velocities

; 7.3.18 Bonds
constraints             = hbonds        ; convert all bonds to constraints
constraint_algorithm    = LINCS         ; LINear Constraint Solver
continuation            = no            ; no = 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
*********************************************************************************
Ubq npt.mdp

; 7.3.2 Preprocessing
define                  = -DPOSRES              ; defines to pass to the
preprocessor

; 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                  = 250000                ; maximum number of steps
to integrate, 0.002 * 250,000 = 500 ps
comm_mode               = Linear                ; remove center of mass
translation
nstcomm                 = 1                     ; [steps] frequency of mass
motion removal

; 7.3.8 Output Control
nstxout                 = 500           ; [steps] freq to write coordinates
to trajectory
nstvout                 = 500           ; [steps] freq to write velocities
to trajectory
nstfout                 = 50000         ; [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
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

; 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

; 7.3.15 Pressure Coupling
pcoupl                  = parrinello-rahman     ; pressure coupling where
box vectors are variable
pcoupltype              = isotropic             ; pressure coupling in
x-y-z directions
refcoord_scaling        = all
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

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

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

; 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

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



More information about the gromacs.org_gmx-users mailing list