[gmx-users] large forces and monstrous water molecules in energy minimization step

[Justin A. Lemkul]

Use a different box size.  I replicated your problem, but your run completes 
successfully with a box set up with editconf -d 1 instead of -d 3.3.

If you set nstxout = 1 during the EM process, you'll see the problematic water 
molecule become unstable.  It looks as if there is a small void in the solvent 
(due to the way genbox tries to stack solvent configurations) and your water 
molecule simply can't find a good orientation within that void.

-Justin

Ehud Schreiber wrote:
> Dear GROMACS users,
> 
>  
> 
> When trying to simulate a pair of interacting proteins in water, I have 
> encountered problems that ultimately resulted in the simulation 
> crashing. I then tried to simplify the system as far as possible while 
> retaining the problem; I now believe the problem (or at least a part of 
> it) lies in the energy minimization step (the first molecular dynamics 
> one). Specifically, the forces encountered during this step are very 
> large, and some water molecules (which are supposed to be rigid) become 
> giant and misshapen.
> 
>  
> 
> In more details:
> 
> 1) I use GROMACS version 4.0.7, single precision, on a server with two 
> Intel x86-64 processors and the redhat 5.4 linux OS.
> 
> 2) I created a PDB file, called NaCl.pdb, with only two “atoms”, 
> actually Na+ and Cl- ions separated by their distance in the lattice of 
> salt:
> 
>  
> 
> HET     NA  A   1       1
> 
> HET     CL  A   1       1
> 
> HETNAM      NA SODIUM ION
> 
> HETNAM      CL CHLORIDE ION
> 
> FORMUL   1   NA    NA 1+
> 
> FORMUL   2   CL    CL 1-
> 
> HETATM    1 NA    NA A   1      -1.41    0.0   0.0    1.00  0.0          NA
> 
> HETATM    2 CL    CL A   1      +1.41    0.0   0.0    1.00  0.0          CL
> 
>  
> 
> 3) I use the tip3p water model:
> 
>  
> 
> pdb2gmx -f NaCl.pdb -water tip3p
> 
>  
> 
> 4) I create the box:
> 
>  
> 
> editconf -f conf.gro -bt dodecahedron -d 3.3 -o box.gro
> 
>  
> 
> 5) I add water using spc216, creating the saltwater.gro file (which 
> seems O.K. by inspection):
> 
>  
> 
> genbox -cp box.gro -cs spc216.gro -p topol.top -o saltwater.gro
> 
>  
> 
> 6) I create the energy minimization parameter file em.mdp:
> 
>  
> 
> ------em.mdp------
> 
> integrator      = steep
> 
> nsteps          = 200
> 
> nstlist         = 10
> 
> rlist           = 1.0
> 
> coulombtype     = pme
> 
> rcoulomb        = 1.0
> 
> vdwtype         = Cut-off
> 
> rvdw            = 1.0
> 
> nstenergy       = 10
> 
> ------------------
> 
>  
> 
> 7) I prepare the em.tpr file for the energy minimization run:
> 
>  
> 
> grompp -f em.mdp -p topol.top -c saltwater.gro -o em.tpr
> 
>  
> 
> 8) I run the energy minimization step:
> 
>  
> 
> mdrun -v -deffnm em
> 
>  
> 
> 9) Looking at the em.log file I see that this step converged to machine 
> precision but did not have maximal force < 10:
> 
>  
> 
> …
> 
> Enabling SPC water optimization for 7561 molecules.
> 
> …
> 
> Max number of connections per atom is 2
> 
> Total number of connections is 30244
> 
> Max number of graph edges per atom is 2
> 
> Total number of graph edges is 30244
> 
> Going to use C-settle (7561 waters)
> 
> wo = 0.333333, wh =0.333333, wohh = 3, rc = 0.075695, ra = 0.0390588
> 
> rb = 0.0195294, rc2 = 0.15139, rone = 1, dHH = 0.15139, dOH = 0.09572
> 
> …
> 
> Stepsize too small, or no change in energy.
> 
> Converged to machine precision,
> 
> but not to the requested precision Fmax < 10
> 
> …
> 
> Steepest Descents converged to machine precision in 36 steps,
> 
> but did not reach the requested Fmax < 10.
> 
> Potential Energy  = -3.4678925e+05
> 
> Maximum force     =  6.4623531e+05 on atom 11052
> 
> Norm of force     =  5.4643726e+03
> 
>  
> 
> 10) Looking at the em.gro file I see one monstrous water molecule (no. 
> 3686); e.g., it has |HW2-OW| = 3.384876 nm, while the normal distance is 
> about 0.1 nm. Its HW2 atom (no. 11054) is close to another water 
> molecule (no. 5849), e.g., 0.047 nm from the latter’s HW2 atom (no. 17543):
> 
>  
> 
> …
> 
> 3686SOL     OW11052   4.348   3.778  -0.629
> 
> 3686SOL    HW111053   5.360   2.601   0.505
> 
> 3686SOL    HW211054   6.518   1.650   0.861  
> 
> …
> 
> 5849SOL     OW17541   6.525   1.698   0.900
> 
> 5849SOL    HW117542   6.606   1.649   0.918
> 
> 5849SOL    HW217543   6.481   1.648   0.832
> 
> …
> 
>  
> 
> 11) During the simulation, several files called stepnnl.pdb were 
> produced for problematic steps, where nn=11,15,19 and l=b,c. For 
> example, the file step19c.pdb indeed shows a problematic water molecule 
> no. 3686, while step19b.pdb does not. Likewise, the earlier step11c.pdb 
> shows a problematic water molecule no. 3266 while step11b.pdb seems 
> proper. The stdout/stderr of mdrun contains warnings like the following:
> 
>  
> 
> …
> 
> t = 0.019 ps: Water molecule starting at atom 11052 can not be settled.
> 
> Check for bad contacts and/or reduce the timestep.
> 
> …
> 
>  
> 
> 12) Reducing the neighbor list update time, i.e., setting nstlist = 1, 
> does not produce any change.
> 
>  
> 
> 13) Trying to use the conjugate gradient integrator instead of steepest 
> descent, i.e., setting integrator = cg, is even worse - the running crashes:
> 
>  
> 
> …
> 
> Step 16, Epot=-1.258771e+35, Fnorm=      nan, Fmax=      inf (atom 14493)
> 
> Segmentation fault
> 
> Exit 139
> 
>  
> 
>  
> 
>  
> 
> So, am I doing something wrong? How can I avoid these problems?
> 
>  
> 
> Thanks,
> 
> Ehud Schreiber.
> 
>  
> 

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

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

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