Measuring the temperature distribution in an MD run

[fgygi]

The following steps describe the extraction of the distribution of temperature for each species during an MD simulations, and the subsequent analysis using the R package through a fit to a gamma distribution.

Assuming that the Qbox output file for the MD simulation is md.r

1) Extract the temperature of a species using qbox_species_temp.sh (provided in the util directory)

Code: Select all

$ qbox_species_temp.sh <species_name> md.r > temp.dat

2) Start R:

Code: Select all

$ R
> library(fitdistrplus)
> t<-read.table("temp.dat")[[1]]

# Fit to a gamma distribution using the maximum likelihood method (mle)
> fit<-fitdist(t,distr="gamma",method="mle")

# print the parameters of the fit
> summary(fit)

# plots of the fit
> plot(fit)

Note: the "fitdistrplus" R package must be installed. As root, use:

Code: Select all

# R
> install.packages('fitdistrplus')

Note: the values of the temperature printed by qbox_species_temp.sh are
computed as:

( vx*vx + vy*vy + vz*vz ) * (2/3)*(1/2)*mass*mass/kB

so that the output is the temperature in K. This allows for comparisons of the
fitted distributions of different species.

The fit parameters are e.g.:

Code: Select all

> summary(fit)
Fitting of the distribution ' gamma ' by maximum likelihood 
Parameters : 
        estimate   Std. Error
shape 1.35513831 1.109630e-02
rate  0.00499023 4.716253e-05
Loglikelihood:  -142092.4   AIC:  284188.8   BIC:  284204.8 
Correlation matrix:
          shape      rate
shape 1.0000000 0.8056883
rate  0.8056883 1.0000000

The shape and rate parameters characterize the gamma distribution.

From the R "help(rgamma)" output:

The Gamma distribution with parameters ‘shape’ = a and ‘scale’ = s
has density

f(x)= 1/(s^a Gamma(a)) x^(a-1) e^-(x/s)

for x >= 0, a > 0 and s > 0. (Here Gamma(a) is the function
implemented by R's ‘gamma()’ and defined in its help. Note that a
= 0 corresponds to the trivial distribution with all mass at point
0.)

The mean and variance are E(X) = a*s and Var(X) = a*s^2.

with scale = 1/rate
In the above example, shape=a=1.35513831 scale=s=1/0.00499023 so that
E(t) = a*s = 1.35513831/0.00499023 = 271.558

i.e. the mean temperature is 271.558 K

The temperature can be computed in R using the fit estimated parameters as

Code: Select all

> fit$estimate[1]/fit$estimate[2]