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

[Ran Friedman]

Hi,

For the dynamics to work you indeed need smaller forces (on the order of
10^3 in GMX units).
Using flexible water molecules I was able to get this for your NaCl model. 
Just add:
define              =  -DFLEXIBLE
To your input.

This should work also to fill in the voids I guess.

Ran

-- 
------------------------------------------------------
Ran Friedman
Postdoctoral Fellow
Computational Structural Biology Group (A. Caflisch)
Department of Biochemistry
University of Zurich
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
Tel. +41-44-6355559
Email: r.friedman at bioc.uzh.ch
Skype: ran.friedman
------------------------------------------------------

Justin A. Lemkul wrote:
>
> 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.