Measuring the temperature distribution in an MD run
Posted: Tue Oct 05, 2021 9:22 pm
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)
2) Start R:
Note: the "fitdistrplus" R package must be installed. As root, use:
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.:
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
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
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)
Code: Select all
# R
> install.packages('fitdistrplus')
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
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]