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

[Justin A. Lemkul]

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

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================