Measuring the temperature distribution in an MD run

Questions and discussions regarding the use of Qbox
Forum rules
You must be a registered user to post in this forum. Registered users may also post new topics if they consider that their subject does not correspond to any topic already present on the forum.
Post Reply
fgygi
Site Admin
Posts: 150
Joined: Tue Jun 17, 2008 7:03 pm

Measuring the temperature distribution in an MD run

Post by 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]
Post Reply