[gmx-users] Box size error during attempted NPT equilibration

[Justin A. Lemkul]

Andrew DeYoung wrote:
> Hi,
> 
> I have tried to run a simulation of 1000 SPC/E water molecules.  I have a
> PDB file containing a 10 by 10 by 10 regular array (cube) of water
> molecules, each separated by approximately 6 Angstroms.  I used pdb2gmx to
> convert the PDB file to GRO and topology, selecting the OPLS-AA force field

Just FYI - it's far easier to create such a topology using a text editor rather 
than pdb2gmx (which is generally not designed to create topologies for multiple 
molecules, rather ones that are chains).  In your case:

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"

[ system ]
water

[ molecules ]
SOL   1000

is all you need.  Maybe that will make your life easier in the future.

> and the SPC/E water model.  I then used editconf to adjust the box size to a
> cube of edge length 60 Angstroms (using the option -box 6.0 6.0 6.0 in

Was the original box incorrectly assigned?  Otherwise, I see no reason to 
intervene with editconf.

> editconf, because the -box option takes arguments of units in nm, I think).

Refer to Chapter 2 for standard units used by all Gromacs programs.  Distances, 
lengths, etc are in nm.

> I chose an edge length of 60 Angstroms because by looking at the GRO file,
> it seems that the largest position coordinate is ~54 Angstroms, so a box
> with edge length 60 Angstroms should, I think, contain all of the 1000 water
> molecules.
> 

You can remove the "I think" part by visualizing the structure in ngmx, VMD (by 
enabling the unit cell representation), etc.

> I then ran energy minimization (steep integrator, emtol = 1000.0 kJ/mol/nm,
> emstep = 0.01 nm, nsteps = 50000), and steepest descents converged to Fmax <
> 1000 in one single step.  The final potential energy was -3333 kJ/mol, the
> final maximum force was 11 kJ/mol/nm, and the final norm of force was 8
> kJ/mol.  It seems a little strange to me that this minimization converged in
> only one step, but I guess that this means that the water molecules were
> placed far enough away from each other to begin with. 
> 

At 6-Angstrom spacing, there's probably not a lot of interaction between the 
molecules.  This seems fine, although the molecules will have to do a LOT of 
rearranging to find an optimal geometry.

> Then I ran an NPT equilibration, and after about 1.5*10^6 steps, the
> simulation crashed and I got the following error message:
> 
> ---
> Program mdrun, VERSION 4.5.4
> Source code file: domdec.c, line: 2861
> 
> Fatal error:
> The X-size of the box (4.000392) times the triclinic skew factor (1.000000)
> is smaller than the number of DD cells (8) times the smallest allowed cell
> size (0.500000)
> ---
> 
> If you have time, do you have any idea about what was the cause of this
> crash?  
> 
> If I understand correctly, the error message seems to be telling me that
> (4.000392*1) < (8*0.5).  Could it be that my box size is fluctuating too
> much during this attempted NPT equilibration, leading to the crash?  Or
> perhaps my value of tau_p (tau_p = 2.0 ps) is too large?
> 

Your box size is definitely fluctuating to much.  You started with a 6-nm cube, 
and now your box is only 4 nm long in the x-dimension.  The system is collapsing 
in on itself.  Larger values of tau_p are generally *more* stable when your 
system is far from equilibrium, because the system does not have to respond as 
quickly to oscillations in the pressure.

> To summarize, I am trying to do a 5 ns NPT equilibration (5*10^6 steps of
> 0.001 ps each) with the md integrator, Nose-Hoover temperature coupling, and
> Berendsen pressure coupling.  Should I next try V-rescale temperature
> coupling instead of Nose-Hoover, since Florian Dommert pointed out to me the
> other day that
> (http://lists.gromacs.org/pipermail/gmx-users/2011-May/061716.html)
> V-rescale tends to be more stable?
> 

Good idea.  My recommendation is to always use either V-rescale or Berendsen 
temperature coupling for initial equilibration, except for the most robust 
systems.  Why not start with a short NVT to help some of the geometry optimize? 
  Right now, your system is trying to re-orient and compress to fill the 
intermolecular voids all at the same time.  That's a lot to ask in one step. 
Doing a bit of NVT (with Berendsen or V-rescale) will allow some of those 
rearrangements to happen more gently.  Applying a barostat to a very 
compressible system (one in which there is a lot of void space) is asking for 
trouble.

> The complete input file that I used for equilibration is as follows:
> ---
> ; npt-eq.mdp
> ; NPT Equilibration (5 ns)
> ; NO POSITION RESTRAINT
> title = OPLS Water NPT equilibration
> 
> integrator = md
> nsteps = 5000000
> dt = 0.001
> 
> nstxout = 250
> nstvout = 250
> nstenergy = 250
> nstlog = 250
> 
> ns_type = grid
> nstlist = 1

You're going to take a major performance hit with nstlist = 1.  A value of 5 or 
10 is much more suitable here.  Setting nstlist = 1 is only necessary for EM.

-Justin

> rlist = 1.0
> rcoulomb = 1.0
> rvdw = 1.0
> 
> coulombtype = PME
> pme_order = 4
> fourierspacing = 0.12
> 
> tcoupl = nose-hoover
> tc-grps = system
> tau_t = 0.1
> ref_t = 298
> 
> pcoupl = berendsen
> pcoupltype = isotropic
> tau_p = 2.0
> ref_p = 1.0
> compressibility = 4.5e-5
> 
> pbc = xyz
> 
> DispCorr = no
> 
> gen_vel = yes
> gen_temp = 298
> gen_seed = 173529
> ---
> 
> Also, the commands I used for equilibration are:
> 
> grompp -f npt-eq.mdp -c em.gro -p water.top -o npt-eq.tpr -po npt-eqout.mdp
> mdrun -v -deffnm npt-eq -cpo npt-eq.cpt
> 
> 
> Thank you very much for your time!  I really appreciate it.
> 
> Andrew DeYoung
> Carnegie Mellon University
> 

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

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

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