[gmx-users] Self-diffusion in high temperature water

[Jose Ignacio Marquez]

Hi,

I came across the following problem: when I calculate the frequency
spectrum for high temperature (523 K) water, the results are inconsistent
with the literature and analytical self-diffusion models.

I am trying to reproduce the calculations by Martí, Guardia and Padro[1] of
the vibrational spectrum of water along the liquid-vapor coexistence line.
To do this, I first equilibrate the configuration for a given temperature
and pressure for 10 o 20 ps with a short (0.1 fs) timestep, and then run
the simulation for another 20 ps with the same timestep, dumping the
velocity to the .trr file every few frames ([2] is the input file). Then,
using g_velacc, I compute the VACF for H, O and the center of mass
velocities.

To compute the frequency spectrum I first compute the VACF for negative
times as VACF(t) = VACF(t), and then calculate the FFT of the simmetrized
VACF.

I have done this for different flexible water models (SPC, SPC/E, TIP4P,
TIP4P/2005) and all seem to reproduce the main features of the frequency
spectrum, with slight variations.

But, for some reason, the shape of the frequency spectra at low energies
differs from what it should be expected for a diffusive regime. Instead of
getting a Lorentzian shape, the limit of the spectrum for omega -> 0 is a
decreasing exponential. At first I thought this could have been caused by
one particular water model, but actually I find the same behavior for all
models, with different rho(omega=0) values, caused by different estimates
of the diffusion coefficient.

At ambient temperature (300 K), the Lorentzian shape is masked by the
presence of intermolecular vibrations (Hydrogen bond bending and
stretching), but at higher temperatures (523 K) diffusion dominate and the
characteristic Lorentzian shape should be seen. To compare, I plotted[3]
the spectrum of center of mass velocity computed at 523K with TIP5P/2005
Flex and the (scaled) spectrum of oxygen velocity, computed by Marti et al.

I suspect the problem might be with tuning the options of GROMACS to ensure
that the long time behavior is correctly represented, but I am fairly new
to GROMACS and MD in general and I don't know how to achieve it. I tried
running for longer times -up to 1 ns-, I tried running on double precision
to ensure that is not a rounding problem, I tried to run on a larger system
(to get a more converged statistical average), I modified g_velacc to save
the VACF with more significant digits... but nothing seems to solve it.

Any help on this will be appreciated.

Thanks in advance,
                Ignacio


[1] http://dx.doi.org/10.1063/1.471932<http://link.aip.org/link/doi/10.1063/1.471932>
[2] http://pastebin.com/z4bC5ZYZ
[3] http://i.imgur.com/53R21.jpg