Hello,
I am trying to calculate the diffusion coefficient of SPCE water from a short MD simulation trajectory. Using the mean square deviation technique, I am able to get the value of D as 2.5991 (1e-5 cm^2/s) which seems to be in the correct ballpark.
However, when I try to derive the value using the velocity autocorrelation function, the value is quite different from that derived using MSD.
I ran a short simulation with the following parameters as the continuation of a 20 ns trajectory.
nsteps = 50000 ;
dt = 0.002
nstxout = 10
nstvout = 10
nstenergy = 10
nstlog = 10
and computed the velocity autocorrelation function using the flag -mol.
Following the steps here: Diffusion Constant - Gromacs, (basically, integrating the area with gmx analyze and dividing by 3)
the calculated value of the diffusion constant comes to around: 5.97 (1e-5 cm^2/s).
I am hoping that someone could point what has changed in how gmx velacc computes the autocorrelation or integrate option of gmx analyze from older versions of gromacs. Or, if there is another way to compute the diffusion coefficient from velocity autocorrelation data.
Why would you want to compute the diffusion constant from the VAC? Using the MSD, which is exactly the integral, is a much better way.
I guess your value is off because the fluctuations in the velocity are very fast, so you need more frequent sampling than every 10 steps. You also don’t state which velocity you are using.
Hello Hess,
Thank you for your response. I simply wanted to check if the VAC would yield a value of diffusion coefficient in the same range as that from MSD. I will stick to the MSD route. However, I am curious about your comment on “which velocity are you using”. I am afraid I don’t understand what you mean by that. When I invoke gmx velacc, I use the -mol flag to signify that I am interested in water molecules as a group. The velocities are taken from the .trr file produced during the simulation. The simulation used the V-rescale thermostat along with Parrinello-Rahman barostat.
The offset could come from the MSD method excluding the ballistic region. With gmx msd it will automatically fit the MSD from 10% to 90%. The first 10 % contain the ballistic region. With the integral of the velocity autocorrelation you always include the ballistic region. But I have to say, with 20 ns (assuming you calculate a 10 ns VACF) this should not play a role. So probably some constants are wrong or the sampling or lenght of the VACF are not practical.
hello~ I’m trying derive the diffusion coefficient using velocity autocorrelation function. I’m puzzled about the output file produced about velacc. Could you please tell me the unit used by C(t). I thought this unit is nm^2/ps^2. Am I right?
