[gmx-users] NVT Equilibration error

[Agnivo Gosai]

Dear Users,

I have a solvated protein-dna structure placed inside a simulation box.
During topology preparation and subsequent steps I did not receive any
warning , except a long bond warning in the pdb2gmx step for the protein.

I then did an energy minimization. The script is shown below:

; For energy minimization
; Parameters describing what to do, when to stop and what to save
integrator  = steep     ; Algorithm (steep = steepest descent minimization)
emtol       = 500.0    ; Stop minimization when the maximum force < 1000.0
kJ/mol/nm
emstep      = 0.01      ; Energy step size
nsteps      = 500000     ; Maximum number of (minimization) steps to perform
energygrps  = System    ; Which energy group(s) to write to disk


; 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.4       ; Cut-off for making neighbor list (short range
forces)
coulombtype = PME       ; Treatment of long range electrostatic interactions
rcoulomb    = 1.4       ; Short-range electrostatic cut-off
rvdw        = 1.4       ; Short-range Van der Waals cut-off
pbc         = xyz       ; Periodic Boundary Conditions

This is the result at the end of the minimization :

Steepest Descents converged to Fmax < 500 in 8938 steps
Potential Energy  = -7.5707590e+06
Maximum force     =  4.4626624e+02 on atom 2581
Norm of force     =  3.2807171e+00

Then I proceeded for a 5 ns NVT equilibration , with the following script :

title       = NVT equilibration
define      = -DPOSRES  ; position restrain the protein and ligand
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 2500000    ; 5 ns
dt          = 0.002     ; 2 fs
; Output control
nstxout     = 500       ; save coordinates every 1 ps
nstvout     = 500       ; save velocities every 1 ps
nstenergy   = 500       ; save energies every 1 ps
nstlog      = 500       ; update log file every 1 ps
nstxtcout   = 1000      ;
; Bond parameters
continuation    = no            ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints
constraints     = all-bonds     ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
ns_type     = grid      ; search neighboring grid cells
nstlist     = 5         ; 10 fs
rlist       = 1.4      ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.4      ; short-range electrostatic cutoff (in nm)
rvdw        = 1.4       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype     = PME       ; Particle Mesh Ewald for long-range
electrostatics
pme_order       = 4         ; cubic interpolation
fourierspacing  = 0.16      ; grid spacing for FFT
; Temperature coupling
tcoupl      = V-rescale                     ; modified Berendsen thermostat
tc-grps     = DNA_Protein Water_and_ions    ; two coupling groups - more
accurate
tau_t       = 0.1     0.1                     ; time constant, in ps
ref_t       = 300     300                     ; reference temperature, one
for each group, in K
; Pressure coupling
pcoupl      = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc         = xyz       ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = yes       ; assign velocities from Maxwell distribution
gen_temp    = 300       ; temperature for Maxwell distribution
gen_seed    = -1        ; generate a random seed
; COM motion removal
; These options remove COM motion of the system
nstcomm         = 10
comm-mode       = Linear
comm-grps       = System

However the simulation crashed and the following error messages were
received :

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

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

-------------------------------------------------------
Program mdrun_mpi, VERSION 4.6.7
Source code file:
/work/gb_lab/agosai/GROMACS/gromacs-4.6.7/src/mdlib/pme.c, line: 851

Fatal error:
1 particles communicated to PME node 11 are more than 2/3 times the cut-off
out of the domain decomposition cell of their charge group in dimension y.
This usually means that your system is not well equilibrated.


Now, what can be going wrong in here ? Do I need to do more energy
minimization , using other algorithms? Or , should I be heating up my
system gradually to 300 K ( from 270 K , for e.g.) and then do a NVT
equilibration at 300 K?

Kindly suggest.


Thanks & Regards
Agnivo Gosai
Grad Student, Iowa State University.