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

[Ehud Schreiber]

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.

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100617/4cda3529/attachment.html>