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]
```