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
========================================
--
gmx-users mailing list [email protected]
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 [email protected].
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
