Statistical thermodynamics

The number of different configurations of positions ${\bf r}^N$ and momenta ${\bf p}^N$ of $N$ particles in a volume $V$ at an absolute temperature $T$ is for all but the simplest systems mindbogglingly large. In the canonical ensemble (i.e. at fixed $N$,$V$, and $T$), the probability $P$ to find the system in some configuration ( ${\bf r}^N, {\bf p}^N$) depends on the total energy of the configuration given by the Hamiltonian ${\cal H}({\bf r}^N, {\bf p}^N)$ according to the Boltzmann distribution:

$$P({\bf r}^N, {\bf p}^N) = Q^{-1} (h^{3N}N!)^{-1} \exp[-\beta {\cal H}({\bf r}^N, {\bf p}^N)]$$ (5)

with $\beta=1/k_BT$. Boltzmann's constant is $k_B=1.38066\cdot10^{-23}$ J/K and Planck's constant is $h=6.62618\cdot10^{-34}$ J$\cdot$s. Here, $Q$ is the partition function given by

$$Q = (h^{3N}N!)^{-1} \int\!\!\int \mathrm{d}{\bf r}^N\!\mathrm{d}{\bf p}^N \exp[-\beta {\cal H}({\bf r}^N, {\bf p}^N)]$$ (6)

and the parenthesized reciprocal term on the right-hand-side of both equantion 2.4 and 2.5 connects the classical distribution to the quantum mechanical distribution. The macroscopic observable $A$ is obtained by taking the ensemble average (denoted by the use of brackets $\left<\cdots \right>$) of the microscopic property $a$,:

$$\left<a\right> = \frac{\int\!\!\int \mathrm{d}{\bf r}^N\!\ma... ...athrm{d}{\bf p}^N \exp[-\beta {\cal H}({\bf r}^N, {\bf p}^N)]}$$ (7)

This horrible expression can numerically be evaluated by sampling the configurations either stochastically, using the Monte Carlo method, or deterministically using molecular dynamics (MD). In the latter case, one makes use of the chaotic (ergodic) nature of the dynamics of the system. That is, during a sufficiently long MD simulation, a representative part of all possible configurations will be sampled, so that the time average of property $a$ equals the naturally weighted ensemble average $\left<a\right>$:

$$\left<a\right> = \lim_{t \to \infty} \frac{1}{t} \int_0^t \mathrm{d}t^\prime a({\bf r}^N, {\bf p}^N)$$ (8)

An advantage of MD over Monte Carlo is that also dynamical properties can be evaluated. Onsager's regression hypothesis states that slow microscopic fluctuations around equilibrium on average decay according to the macroscopic laws.[3] For example, transport properties such as the diffusion coefficient of some species $i$, $D_i$, which is the macroscopic proportionality constant between the flux $J_i$ of the diffusing species $i$ and the gradient of its concentration $c_i$ as given by Fick's law,

$$J_i = -D_i \nabla c_i$$ (9)

can be calculated from an equilibrium MD trajectory (i.e. without inducing a gradient in the concentration) from the microscopic fluctuations in the velocity of the species, when averaged over long enough times $\tau$:

$$D_i = \int_0^\infty \mathrm{d}\tau \left< v_i(\tau)v_i(0)\right>$$ (10)

Simulation of a chemical reaction and therefore the direct estimation of the rate constant $k$ is in practice not possible using MD or Monte Carlo. The problem is, that for the system to move from the stable reactant state to the stable product state, it has to cross some transition state configuration associated with the relatively high activation energy, $\Delta E^\ddagger \gg k_BT$, of equation 2.2, which has a very low probability of being populated during a simulation, as follows from the exponential dependence in equation 2.4. This dynamical bottleneck makes the chemical reaction a rare event on the time scale of the thermal motions. Fortunately, there are a number of techniques available to circumvent this problem.