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

[Jose Ignacio Marquez]

Sorry for the bump, but, please, does anybody know something of any study
of the self-diffusion frequency spectrum of water done with GROMACS?
I am not sure if it is something I am doing wrong, or if that is simple not
possible to model with the program.

Thanks,
         Ignacio

On Tue, Oct 30, 2012 at 6:05 PM, Jose Ignacio Marquez <
nachomarquez at gmail.com> wrote:

> 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
>
>