[gmx-users] problem with water models in Gromacs 5.1.2

[Irem Altan]

Hi,

It turns out that my problems with TIP5P were more general than I thought. I recently switched from version 4.6.5, with which I could run simulations using AMBER99SB + (SPC, SPC/E, TIP4P, TIP3P). I decided to do the following test:

1- I created a .pdb file with a single calcium ion, within a cubic unit cell of side 30 A:

(test.pdb)
CRYST1   30.000   30.000   30.000  90.00  90.00  90.00 P 1           1
HETATM    1 CA    CA A 1        32.702 -10.496  -5.429  1.00 12.22          CA

2- I then run pdb2gmx and select amber99sb and SPC (I tried it with TIP4P as well):

gmx pdb2gmx -f test.pdb -o box.gro -p topol.top

3- I run solvate:

gmx solvate -cp box.gro -cs spc216 -p topol.top -o box_h2o.gro

4- then I substitute some waters with sodium and chloride ions:

gmx grompp -f ions.mdp -c box_h2o.gro -o ions.tpr -p topol.top
gmx genion -s ions.tpr -p topol.top -o ions.gro -neutral -conc 0.05

5- Energy minimization (this is where things go horribly wrong):

gmx grompp -f minim.mdp -p topol.top -c ions.gro -o em.tpr
gmx mdrun -v -deffnm em

Then the following error appears:

Fatal error:

step 11: Water molecule starting at atom 86 can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.

I’m guessing there is something that changed from version 4.6.5 to 5.1.2 that I can’t pinpoint. What could I be overlooking? (The contents of my .mdp files are below)

Best,
Irem

—————

——————
ions.mdp
——————

 ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep         ; Algorithm (steep = steepest descent minimization)
emtol           = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep          = 0.01          ; Energy step size
nsteps          = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1             ; Frequency to update the neighbor list and long range forces
ns_type         = grid          ; Method to determine neighbor list (simple, grid)
rlist           = 1.1           ; Cut-off for making neighbor list (short range forces)
coulombtype     = PME           ; Treatment of long range electrostatic interactions
rcoulomb        = 1.1           ; Short-range electrostatic cut-off
rvdw            = 1.1           ; Short-range Van der Waals cut-off
pbc             = xyz           ; Periodic Boundary Conditions (yes/no)

——————
minim.mdp
——————

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep         ; Algorithm (steep = steepest descent minimization)
emtol           = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep          = 0.01          ; Energy step size
nsteps          = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1             ; Frequency to update the neighbor list and long range forces
ns_type         = grid          ; Method to determine neighbor list (simple, grid)
rlist           = 1.0           ; Cut-off for making neighbor list (short range forces)
coulombtype     = PME           ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0           ; Short-range electrostatic cut-off
rvdw            = 1.0           ; Short-range Van der Waals cut-off
pbc             = xyz           ; Periodic Boundary Conditions (yes/no)