[gmx-users] regarding range checking error and constraint errors in Lincs algorithm

[shyamala iyer]

Hi gromacs users,

I apologize in advance for the long mail, but I have been having some
trouble with my gromacs simulations lately.

My system consists of two popypeptide chains and a small ligand (or
drug like compound). I obtain the required gromacs topologies for the
ligand/drug from PRODRG dundee server. The steps I use for a MD
simulation using gromacs is briefly as follows:
1) get ligand file with required polar hydrogens and gromacs ligand
.itp file from prodrg server
2) run pdb2gmx using the gmx force field and spce water
3) correctly update the protein topology file to include ligand itp
file header, number of ligand molecules and also make sure the output
file from pdb2gmx contains the necessary ligand coordinates
3) set up box around this system, fill box with solvent (water),
neutralize system using genion
4) set up files to run energy minimization on system, run energy minimization
5) set up files to run position restrained dynamics and run position
restrained dynamics
6) setup up files and run Molecular dynamics on the system for ~ 10 ps
The above set up worked with one of the protease-ligand complex I was
testing, but had been failing for another viral protease.

The errors vary depending on the run, and in some simulations, I get
an error during the Energy minimization step: (section of em log file
pasted below)
------------------------------------------------------------


Enabling SPC water optimization for 41607 molecules.

Will do PME sum in reciprocal space.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++

U. Essman, L. Perela, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen

A smooth particle mesh Ewald method

J. Chem. Phys. 103 (1995) pp. 8577-8592

-------- -------- --- Thank You --- -------- --------

Removing pbc first time

Done rmpbc

Initiating Steepest Descents

Center of mass motion removal mode is Linear

We have the following groups for center of mass motion removal:

0: rest, initial mass: 129215

Started Steepest Descents on node 0 Wed Dec 19 10:33:25 2007

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++

S. Miyamoto and P. A. Kollman

SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithms for
Rigid Water Models

J. Comp. Chem. 13 (1992) pp. 952-962

-------- -------- --- Thank You --- -------- --------

Steepest Descents:

Tolerance (Fmax) = 1.00000e+02

Number of steps = 100

Going to use C-settle (41607 waters)

wo = 0.333333, wh =0.333333, wohh = 3, rc = 0.08165, ra = 0.0384897

rb = 0.0192449, rc2 = 0.1633, rone = 1, dHH = 0.1633, dOH = 0.1

Grid: 15 x 15 x 15 cells

Configuring nonbonded kernels...

Testing x86_64 SSE support... present.

Step Time Lambda

0 0.00000 0.00000

-------------------------------------------------------

Program mdrun, VERSION 3.3.1

Source code file: nsgrid.c, line: 226

Range checking error:

Explanation: During neighborsearching, we assign each particle to a
grid based on its coordinates. If your system contains collisions or
parameter

errors that give particles very high velocities you might end up with
some coordinates being +-Infinity or NaN (not-a-number). Obviously, we
cannot put these on a grid, so this is usually where we detect those
errors.

Make sure your system is properly energy-minimized and that the
potential energy seems reasonable before trying again.

Variable ci has value -2147483648. It should have been within [ 0 .. 3375 ]

Please report this to the mailing list (gmx-users at gromacs.org)

-------------------------------------------------------

At the end of this mail I have also pasted the em.mdp file I use to
set up energy minimizations. I would really appreciate it if someone
could help me address this Range checking error issue.

Another type of error that I recently encountered occurs during
position restrained dynamics simulations. The error is a constraint
error in Lincs algorithm at a particular time step and usually is
accompanied by a line that is similar to, (with time and atom number
varying )
-----------------------------------------------------------------------------------------------------------------------------------------
t= 0.02ps Water molecule starting at atom 71500 can not be settled.
Check for bad contacts and/or reduce the timestep.Wrote pdb files with
previous and current coordinates
-------------------------------------------------------------------------------------------------------------------------------------------
I have triend reducing the time step and I still get the same error.
I have also pasted the pr.mdp file I use to set up and run position
restrained dynamics.
I would like some help in solving this lincs error issue.

The original protease pdb file that I used for the above simulations
has a small segment missing from chain A. This loop segment was highly
disorderd in the cystal structure hence, coordinates for this 7
residue segment is missing in the pdb file. This missing segment is no
where near where the ligand is bound on the protease. Does a missing
segment cause problems during siumlations?

Another general question on the genion function:
Is there a way I could not have genion be interactive. My systems
usually have a protein, a ligand and water, and I always choose group
13 (SOL) as the selection after running genion. I am asking this
because there are times when I want to run the same protein with
several different ligands and I would like a script to run through all
gromacs steps and not have a particular step be interactive.


Thanks for reading a long mail.
Shyamala
------------------------------------------------------
em.mdp
title               =  fws
cpp                 =  /usr/bin/cpp ; location of cpp
define              =  -DFLEX_SPC
constraints         =  none
integrator          =  steep
dt                  =  0.001    ; ps !
nsteps              =  100 ; 100 steps of energy min
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
coulombtype         =  PME
rcoulomb            =  1.0
vdwtype             =  cut-off
rvdw                =  1.4
table-extension	    =  3.0
fourierspacing		=  0.12
fourier_nx		=  0
fourier_ny		=  0
fourier_nz		=  0
pme_order		=  6
ewald_rtol		=  1e-5
optimize_fft		=  yes
;
;       Energy minimizing stuff
;
emtol               =  100.0
emstep              =  0.01
-------------------------------------------------------
pr.mdp
title               =  fws
cpp                 =  /usr/bin/cpp
define              =  -DPOSRES
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002	; ps !
nsteps              =  5000	; total 10.0 ps.
nstcomm             =  1
nstxout             =  250
nstvout             =  1000
nstfout             =  0
nstlog              =  10
nstenergy           =  10
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
coulombtype         =  PME
rcoulomb            =  1.0
vdwtype             =  cut-off
rvdw                =  1.4
fourierspacing		=  0.12
fourier_nx		=  0
fourier_ny		=  0
fourier_nz		=  0
pme_order		=  6
ewald_rtol		=  1e-5
optimize_fft		=  yes
; Berendsen temperature coupling is on
Tcoupl              =  berendsen
tau_t               =  0.1	0.1
tc_grps		    =  protein	non-protein
ref_t               =  300	300
; Pressure coupling is on
Pcoupl              =  berendsen
pcoupltype          =  isotropic
tau_p               =  1.0
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is on at 300 K.
gen_vel             =  yes
gen_temp            =  300.0
gen_seed            =  173529