[gmx-users] g_velacc again

[Nuno R. L. Ferreira]

Dear all gmx users

Few days ago I asked about some hints on how to get diffusion coefficients
from VACF (Cvv(t)).

Read some stuff about the subject (Allen & Tildesley, tkx Erik) and in other
net sources.
As Anton said :

<   the integral of your VACF is exactly the same as the slope of the
positional MSD (or something similar).

If I'm not wrong,  "slope of the positional MSD" = Dt (self-diffusion
coefficient) = 1/3 * (integral of VACF) , speaking in 3D.

I've done a 100 ps MD, 2 fs of step and saving velocities every step (tkx
Anton), in a methanol box (400 molecules, tkx Christoph Freudenberger), and
tryed to compute the diffusion coefficient from MSD (g_msd) and also from
VACF (g_velacc). The Cvv(t) vs time graph shows an expected trend for this
kind of correlations (it decreases fast to some negative value, and then a
Cvv tail aproaches to zero)

g_msd gives a good value ( 2.7 +- .2 10-5 cm2 s-1), compared with the
experimental one ( 2.4 10-5 cm2 s-1).
But I "cannot manage" to obtain the value from g_velacc calculation. In
fact, I don't understand very well the output from g_velacc.
I've employed this script with the following flags:

g_velacc -f traj.trr -s traj.tpr -n

being the index file 10 methanol molecules to speed up things.
Here's the output:
COR: Correlation time (plain integral from 0.000 to 99.999 ps) = 0.01533 ps

What is this correlation time?
With -fitfn aexp , I obtain the same information as above plus:
COR: Relaxation times are computed as fit to an exponential:

COR: y = a2 exp(-x/a1)

COR: Fit to correlation function from 0.000 ps to 100.001 ps, results in a

COR: Fit from            Integral         Tail Value          Sum (ps)
a1 (ps)            a2 ()

COR: 0.0000e+00   0.0000e+00   6.0026e-02   6.0026e-02   5.1511e-02
1.1653e+00

I was expecting to get an integral value from here, but it's zero. I'm also
aware that the tail in Cvv(f) function as to be treated with special care
( Solvent diffusion outside macromolecular surfaces, Eric Lindahl & Olle
Edholm, Phys Rev. E ,V57, N1, Y1998, P791).
I used this function (aexp) cause I read somewhere (actually in a NAMD
tutorial, last summer course) that we can fit the Cvv(t) to a first order
exponential, like this :

Cvv(t) = Dt / tau * exp (-t / tau)
This example was applyed to ubiquitin in water, to calculate the Dt.

Can anibody give me some help on this? My aim is to obtain the diffusion
coefficient for a peptide. Since that for calculating MSD's of proteins with
need longer simulations, and the statisticall error with just one molecule
is "big" (MSD vs time does not show a good straight line), I was thinking in
VACF vs time. But till now I did not understood how to calculate Dt from the
information of g_velacc calculation. Where is the integral of VACF? Am I
missing something?

Many thanks,
Nuno

P.S. I apologise if my explanation is a bit confuse.














#######
Nuno Ricardo Santos Loureiro da Silva Ferreira
Grupo de Química Biológica
Departamento de Química
Faculdade de Ciências e Tecnologia
Universidade de Coimbra
3004-535 Coimbra
Portugal

Phone: +351 239 852080
Fax: +351 239 827703
www.biolchem.qui.uc.pt
####