[gmx-users] Velocity autocorrelation function of water using g_velacc

[David van der Spoel]

On 2011-05-23 20.52, Justin A. Lemkul wrote:
>
>
> Andrew DeYoung wrote:
>> Hi,
>> I have run a simulation of 254 SPC/E water molecules in a cubic box of
>> edge
>> length 19.660 Angstrom. I have used the gromacs utility g_velacc to
>> calculate the velocity autocorrelation function for the motion of only
>> the
>> oxygen atoms. When I plot the velocity autocorrelation function versus
>> time, my result is a curve with more than one local minimum.
>> My question is, should I expect having more than one local minimum? I am
>> not sure that my result is reasonable.
>>
>> For example, the Democritus molecular dynamics website gives one
>> example of
>> a velocity autocorrelation function, although they do not specifically
>> say
>> what system (water or something else?) or what state (liquid or gas?) was
>> used:
>> http://www.cse.scitech.ac.uk/ccg/software/Democritus/Theory/vaf.html .
>> Their plot shows one "strong" local minimum and one "weaker" one; both of
>> these local minima are in the negative regime. I am new to this field,
>> but
>> I believe that the "strong" local minimum corresponds conceptually to the
>> recoil of the atoms when they collide.
>> If you have time, could you please look at my result at
>> http://www.andrew.cmu.edu/user/adeyoung/velacc_may23.pdf and see if you
>> think it is a qualitatively reasonable result? My main concern is the
>> first
>> local minimum, which I have labeled with an arrow on the plot. This
>> minimum
>> is in the _positive_ regime, which I believe (but I am not sure) cannot
>> correspond to the recoil of atoms when they collide.
>> In my pdf file I have displayed the parameters I have used for
>> equilibration
>> and dynamics. To summarize, I have used the leap frog algorithm with 2 ns
>> equilibration in 1 fs steps, and 5 ns dynamics in 1 fs steps. For
>> dynamics,
>> I have saved configurations every 2 fs. Although my result from g_velacc
>> has velocity autocorrelation frunction data from 0 to 1000 ps, the
>> function
>> decays rapidly to zero (as it should, I think, because the velocity
>> rapidly
>> decorrelates since the atoms/molecules "forget" their initial
>> velocities),
>> so I have plotted only from 0 to 1 ps. I am running Gromacs version
>> 4.5.4.
>>
>>
>
> I don't know if this explains your result or not, but you do realize
> that if you have "gen_vel = yes" at the beginning of your simulation for
> data collection, it negates any value in previous equilibration, right?
> You've destroyed the ensemble you created before. You should pass the
> checkpoint file from the equilibration run to grompp -t and set "gen_vel
> = no" for your simulation. I'd start by doing that and seeing if it
> affects the result, because at this point, it may just be that you
> haven't properly re-equilibrated, or at the very least, the original
> frames are messing with the VACF calculation.
>
> -Justin
>
Or you can just omit the first 100 ps from the analysis.
There can be multiple minima, not a problem. If you integrate the VACF 
you should get the diffusion constant (3D), check manual and compare to 
g_msd.

-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se