[gmx-users] cool down problem

[Riccardo Concu]

Dear all, I'm having some problems during the cool-down process after
the NVT. 
I'm running the NVT process at a T of 500K for 3ns and everything is
going well but when I cool down the system blows up. I am running the
cool down for 100ps, I tried to run the process for longer time (500 ps)
and decrease nstlist option up to 5 but it seems it doesn't work. The EM
seems to be fine with very low forces:
Potential Energy  = -3.1736706e+05
Maximum force     =  2.6573439e+01 on atom 105
Norm of force     =  3.1027287e-01

The NVT is compiled with the following mdp:
define      = -DPOSRES  ; position restrain the protein and ligand
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 1500000     ; 2 * 50000 = 100 ps
dt          = 0.002     ; 2 fs
; Output control
nstxout     = 10000       ; save coordinates every 0.2 ps
nstvout     = 10000       ; save velocities every 0.2 ps
nstenergy   = 10000       ; save energies every 0.2 ps
nstlog      = 10000       ; update log file every 0.2 ps
energygrps  = NAP NA WAT MET SI3 SI- I LI+ LI-
; 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
vdwtype = cutoff
ns_type     = grid      ; search neighboring grid cells
nstlist     = 20         ; 10 fs
rlist       = 1.0       ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw        = 1.0       ; 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     = NAP NA WAT MET SI3 SI- I LI+ LI-          ; two coupling
groups - more accurate
tau_t       = 0.1   0.1  0.1  0.1  0.1  0.1     0.1 0.1
0.1                     ; time constant, in ps
ref_t       = 500 500 500 500 500 500 500 500 500          ; 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    = 298       ; temperature for Maxwell distribution
gen_seed    = -1        ; generate a random seed

The cooldown mdp is the following:

integrator               = md
tinit                    = 0
dt                       = 0.001
nsteps                   = 100000; 100 ps
simulation_part          = 1
init_step                = 0
comm-mode                = Linear
nstcomm                  = 1
comm-grps                = 
bd-fric                  = 0
ld-seed                  = 1993

rtpi                     = 0.05
nstxout                  = 0
nstvout                  = 0
nstfout                  = 0
nstlog                   = 500
nstenergy                = 100
nstxtcout                = 500
xtc-precision            = 1000
xtc-grps                 = 
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
periodic_molecules       = no
rlist                    = 0.9
coulombtype              = PME
rcoulomb-switch          = 0
rcoulomb                 = 0.9
epsilon_r                = 1
epsilon_rf               = 1
vdw-type                 = Cut-off
rvdw-switch              = 0
rvdw                     = 1.0
DispCorr                 = No
table-extension          = 1
fourierspacing           = 0.12
fourier_nx               = 0
fourier_ny               = 0
fourier_nz               = 0
pme_order                = 6
ewald_rtol               = 1e-5
ewald_geometry           = 3d
epsilon_surface          = 0
optimize_fft             = yes

gb_algorithm             = Still
nstgbradii               = 1
rgbradii                 = 2
gb_epsilon_solvent       = 80
gb_saltconc              = 0
OBC(II)
gb_obc_alpha             = 1
gb_obc_beta              = 0.8
gb_obc_gamma             = 4.85
sa_surface_tension       = 2.092

Tcoupl                   = v-rescale
tc-grps                  = system
tau_t                    = 0.1
ref_t                    = 1000
Pcoupl                   = Parrinello-Rahman
Pcoupltype               = Isotropic
tau_p                    = 2
compressibility          = 4.5e-5
ref_p                    = 1
refcoord_scaling         = No
andersen_seed            = 815131
annealing                = single
annealing_npoints        = 2
annealing_time           = 0 80
annealing_temp           = 500 298
gen_vel                  = no
gen_temp                 = 500
gen_seed                 = 173529
constraints              = all-bonds
constraint-algorithm     = Lincs
continuation             = no
iterations
Shake-SOR                = no
shake-tol                = 0.0001
lincs-order              = 4
lincs-iter               = 1
lincs-warnangle          = 30
morse                    = no

excluded
energygrp_excl           = 

nwall                    = 0
wall_type                = 9-3
wall_r_linpot            = -1
wall_atomtype            = 
wall_density             = 
wall_ewald_zfac          = 3

pull                     = no

disre                    = No
disre-weighting          = Conservative
disre-mixed              = no
disre-fc                 = 1000
disre-tau                = 0
nstdisreout              = 100
orire                    = no
orire-fc                 = 0
orire-tau                = 0
orire-fitgrp             = 
nstorireout              = 100
dihre                    = no
dihre-fc                 = 1000

free-energy              = no
init-lambda              = 0
delta-lambda             = 0
sc-alpha                 = 0
sc-power                 = 0
sc-sigma                 = 0.3
couple-moltype           = 
couple-lambda0           = vdw-q
couple-lambda1           = vdw-q
couple-intramol          = no
acc-grps                 = 
accelerate               = 
freezegrps               = 
freezedim                = 
cos-acceleration         = 0
deform                   = 
E-x                      = 
E-xt                     = 
E-y                      = 
E-yt                     = 
E-z                      = 
E-zt                     = 
user1-grps               = 
user2-grps               = 
userint1                 = 0
userint2                 = 0
userint3                 = 0
userint4                 = 0
userreal1                = 0
userreal2                = 0
userreal3                = 0
userreal4                = 0



Thank you all in advance,
Riccardo