I've been updating my GROMACS Suite from 4.0.7. to 4.5.4 and suddenly have problems with my bulk water (added via genbox) near the periodic boundaries.
I try to minimize my solvated protein (150 mM NaCl, TIP4P Water model, OPLS/AA FF) via steepest descent (Coulomb: Particle Mesh Ewald) with pbc = xyz (worked perfectly for me in 4.0.7). Minimization crashes after several steps and GROMACS tells me that a charge group traveled too far in the y-direction. Visual inspection tells me that the affected molecule is bulk water near the periodic boundary and has no steric clashes or something. I try to delete the water molecule, then another water molecule near the pbc is affected and solvated the protein again several times - no effect. Did I mess something up in my .mdp file? Are there new arguments I have to add in 4.5? Can anybody help?
Why have you chosen to get your solvated box 150 mM NaCl? (instead to add just few ion particles in order to neutralize your system).
I mean, I don't know at what extent the ionic strength is well reproduced and, in the meanwhile, what computational artifacts you may introduce.
As for artifacts, there are also numerous examples for the importance of explicit ions in the literature - In my experience a solute that is rather physiological gives quite reasonable results, at least for my protein. I've also checked the direct environment of my water molecule for ions that could cause the observed problems ... they're quite far away so I don't think they're responsible for the crash
Dear Mr. mann
Try using my mdp file. I used it in Gromacs 4.5 without any problem.
best,
HHA
@Daniel - can you post the .mdp file you are using?
@Darrio - ionic strength in simulations can be very critical, though some force fields have known defects. Simply neutralizing the system with counterions is common, but I have seen firsthand that increasing ionic strength can change a proteins' behavior (in a physiologically relevant way!) and for nucleic acids and membranes, the effects are even more pronounced.
@Hamed - that .mdp file would be rather unconventional for OPLS. Typically, cutoffs of 1.0 nm are used, though there is some dispute about what is "best" for OPLS.
Hey Justin,
here's the used .mdp file for energy minimization:
title = xxx
cpp = /lib/cpp
define = -DFLEXIBLE
integrator = steep
nsteps = 75000
nstxout = 2000
nstvout = 2000
nstfout = 2000
nstlog = 2000
nstenergy = 2000
nstxtcout = 2000
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 0.9
coulombtype = PME
rcoulomb = 0.9
vdwtype = cut-off
rvdw = 1.4
fourierspacing = 0.12
PME_order = 4
ewald_rtol = 1e-5
emtol = 0.001
emstep = 0.01
nstcgsteep = 10
energygrps = SOL NA CL
For minimization, nstlist should always be set to 1. You're probably missing interactions with the updates, leading to instability. The fact that it "worked" with a previous version is probably a happy coincidence, but nstlist would be my #1 suspect at this point.
- Badges
- Science topic