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
========================================
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
