[gmx-users] Exploding box after 55.5ns of simulation

[JW Feng]

Dear GROMACS users,

I am running a NPT simulation of a 28 residue protein in a truncated
dodecahedron box of TIP4P water at 300K and 1bar.  The force field is
OPLSAA 2001.  The system ran fine for 55.5 ns and then it crashed due
to an exploding box.  As the box was exploding, the potential energy
increased from -209 kJ to -204 kJ, the box volume went from 150 to
223, the pressure shot up to 897 bar from an average of 1bar, the
temperature increased from 300K to 317 K, and the total energy
increased from -171 kJ to -164 kJ.

>From reading the mailing archive, perhaps changing tau_p from 1.0 to
10.0 or 20.0 may solve the problem.  Can anyone else please provide
suggestions?  Others exploding box problems documented in the mailing
archive seem to occur early in the simulation and not after 55.5 ns.
It troubles me that the system ran fine for 55.5 ns before crashing.
Starting at time point 41.2 ns, "Correcting invalid box:" messages
start to appear but that's harmless according to various posting in
the mailing list.

This is how I generated my box:
editconf -f Pro8.gro -o Pro8.gro -d 1.2 -bt dodecahedron
genbox -cp Pro8.gro -cs tip4p.gro -o Pro8_b4em.gro -p Pro8.top

I am using gromacs 3.3.1  I have provided a snippet of the log that
warns of exploding box and my mdp file.
Correcting invalid box:
old box (3x3):
   old box[    0]={ 7.11699e+00,  0.00000e+00,  0.00000e+00}
   old box[    1]={ 0.00000e+00,  7.11699e+00,  0.00000e+00}
   old box[    2]={ 3.56227e+00, -3.55472e+00,  4.98818e+00}
new box (3x3):
   new box[    0]={ 7.11699e+00,  0.00000e+00,  0.00000e+00}
   new box[    1]={ 0.00000e+00,  7.11699e+00,  0.00000e+00}
   new box[    2]={-3.55472e+00, -3.55472e+00,  4.98818e+00}
Grid: 14 x 14 x 10 cells
WARNING: your box is exploding! (ncells = 1960)
           Step           Time         Lambda
       27750000    55500.00391        0.00000

Below is my mdp file
title           = Pro
cpp             = /lib/cpp
integrator      = md
dt              = 0.002
nsteps          = 150000000 ; 300ns
nstxout         = 10000
nstvout         = 10000
nstlog          = 10000
nstenergy       = 10000
nstxtcout       = 1000
nstcomm         = 1
comm_grps       = Protein
xtc_grps        = protein
energygrps      = protein CL- SOL
ns_type         = grid
coulombtype     = pme
nstlist         = 10
rlist           = 1.0
rvdw            = 1.0
optimize_fft    = yes
tcoupl          = berendsen
tc-grps         = protein SOL CL-
tau_t           = 0.1 0.1  0.1
ref_t           = 300 300 300
Pcoupl          = Parrinello-Rahman
tau_p           =  1.0
compressibility =  4.5e-5
ref_p           =  1.0
gen_vel         = yes
gen_temp        = 300
gen_seed        = -1
constraints     = hbonds

Thank you,

JW

-- 
JW Feng
Ph.D. Candidate
Center for Computational Biology
Washington University in St. Louis