[gmx-users] Atoms get frozen with "Nose Hoover thermostat with Parrinello-Rahman barostat" for a system of an ion of charge +2 in flexible water molecules

[Surya Prakash Tiwari]

Dear Gromacs users,

I am having a very strange problem with "Nose Hoover thermostat with
Parrinello-Rahman barostat" NPT simulations for a system of an ion of
charge +2 in flexible water molecules. Flexible water is taken from J.
Chem. Phys. 124, 024503 (2006); http://dx.doi.org/10.1063/1.2136877.
The ion is [UO2]2+ with charge on U=2.5 and each O has -0.25.

Atoms very soon get frozen after the simulation starts and they remain
frozen, they do not move at all. Only simulation box oscillates, which
causes atoms to oscillate little bit but they do not move at all.

Pressure, temperature etc. all get messed up.
Here is some output from g_energy:
"
Statistics over 400001 steps [ 0.0000 through 400.0000 ps ], 3 data sets
All statistics are over 40001 points

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Potential                  -23648.5         20    2891.04   -117.711  (kJ/mol)
Temperature                 306.684        0.1    108.192  -0.566776  (K)
Pressure                    538.625        4.4    89041.8     79.925  (bar)
"

I am using this "Nose Hoover thermostat with Parrinello-Rahman
barostat" after a well equilibration configuration obtained with
Berendsen's thermostat and barostat.

This problem is happening only when Parrinello-Rahman barostat used
with Nose Hoover thermostat.

I also noted the following observations:
If I change either one of the thermostat or barostat with Berendsen,
frozen problem disappears.
If I change flexible water to rigid water, problem disappears.
If my system has no ion, no problem happens.
If I play with cut offs, rlist, then also sometimes problem disappears.
The problem is independent on the size of box and number of water molecules.
As expected, Diffusion constant is very close to zero.


Here is my input MDP file:
integrator  = md
dt          = 0.001
tinit       = 0
nsteps      = 400000
nstcomm     = 10
pbc         = xyz


;constraint_algorithm    = lincs
;constraints             = all-bonds
continuation            = yes       ; continuing from equilibrated NPT

nstxout                  = 200
nstvout                  = 50000
nstfout                  = 0
nstlog                   = 5000
nstenergy                = 2000
nstxtcout                = 0
xtc-precision            = 1000


nstlist                  = 10
ns_type                  = grid

rlist                    = 1.1
coulombtype              = PME_switch
rcoulomb_switch          = 0.6
rcoulomb                 = 0.9
vdw-type                 = switch
rvdw-switch              = 0.6
rvdw                     = 0.9

DispCorr    = EnerPres


fourierspacing  = 0.12
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
pme_order       = 6
ewald_rtol      = 1e-6
optimize_fft    = yes

tcoupl      = nose-hoover
tc-grps     = System
tau_t       = 2
ref_t       = 298.15
pcoupl              = Parrinello-Rahman
pcoupltype          = isotropic
tau_p               = 2.0           ; time constant, in ps
ref_p               = 1.0           ; reference pressure (in bar)
compressibility     = 4.5e-5        ; isothermal compressibility, bar^-1

gen_vel                  = no


Any help will be highly appreciated. Please let me know if someone
needs any additional information to solve this issue.

Thanks,
Surya Prakash Tiwari