3.3. Statistical Thermodynamics

---

download:pdf

---

The purpose of this chapter, which is largely inspired by a textbook by Ungerer, Tavitian & Boutin [1] is to help MedeA GIBBS and MedeA LAMMPS users to understand the link between forcefield based molecular simulation methods (Monte Carlo, molecular dynamics) and classical thermodynamics.

In contrast to quantum chemistry, the detailed electronic structure of each molecule is not considered in molecular simulation with classical forcefields. Indeed, such detailed treatment is generally not necessary for thermodynamic properties, as the related molecular interactions do not involve major modifications of electronic structure. In compensation, molecular simulation can consider larger systems, and this often makes it possible to derive macroscopic properties (i.e. that can be compared with measured quantities).

For systems at equilibrium, this connection is provided by the well-established framework of statistical thermodynamics [2] which we recall briefly in section Statistical ensembles and partition functions As the distribution of energy among molecules plays a central role in statistical mechanics, we will then review the way potential energy is modeled in forcefields, see section Forcefields for Materials Simulations and section Monte Carlo Methods. Monte Carlo Methods introduces the various types of Monte Carlo algorithms that are implemented in MedeA GIBBS, while more practical aspects of simulation and numerical methods are addressed in Morphology.

3.3.1. Statistical Ensembles and Partition Functions

The link between microscopic and macroscopic properties is not intuitive, because fluctuations of properties are significant in microscopic systems. It is therefore required to collect multiple snapshots of a microscopic system to compute a meaningful average property. In statistical thermodynamics, this collection of snapshots is called a statistical ensemble.

A statistical ensemble is a collection of various states of the system, which differ in positions and velocities of the component particles. The space of all possible system states, which is of dimension 6N for N particles, is called the phase space.

There are several different types of statistical ensemble, depending on the type of real system we are investigating and the conditions in which it is placed.

3.3.1.1. Canonical or NVT Ensemble

A system at imposed volume V and temperature T, is represented by a statistical ensemble, which is called the canonical ensemble or NVT ensemble. The number of molecules N and the global volume V are identical for all system states belonging to the ensemble, but they differ in total energy E which is a fluctuating variable in this ensemble. Each state j of the canonical ensemble occurs with a probability proportional to \(exp (-E_{j}/k_{B}T)\) where \(k_{B}\) is the Boltzmann constant (\(k_{B} = R/N_{A} , R=8.314J/mol/K\) the gas constant and \(N_{A}=6.0222 \cdot 10^{23}\) the Avogadro number) and \(E_{j}\) is the total energy (kinetic + potential) of the system in state j. Here we use the symbols \(E\) for total energy (kinetic and potential) and \(U\) for the potential energy.

The Boltzmann factor \(exp (-E_{j}/k_{B}T)\) expresses that low energy states are favored compared with high energy states. Increasing temperature broadens the energy distribution in the ensemble with the consequence that the average energy is increased. A general notation used in statistical mechanics is:

\[\beta = \dfrac{1}{k_{B}T}\]

Thus, the probability of a given state j existing in the phase space is given by:

\[P_{j} = \dfrac{exp(\beta E_{j})}{Q_{NVT}}\]

Where

\[Q_{NVT} = \dfrac{1}{N!}\sum_{r_{i}}\sum_{p_{i}} exp \left( -\beta E (r_{i},p_{i}) \right)\]

The expression Q_NVT, known as the partition function, is the sum of the Boltzmann factors for all possible different states in the phase space. The factor 1/N! originates from the possible combinations of N identical particles which correspond to the same state. Consistently with quantum mechanics, the finite summation may be transformed in a continuous integral:

\[Q_{NVT} = \dfrac{1}{h^{3N}N!} \int_{r_i} \int_{p_{i}} exp \left( -\beta E(r_{i},p_{i}) \right) dr_{i} dp_{i}\]

where the factor \(1/h^{3N}\) may be understood as the volume of an individual quantum state in the phase space (involving the Planck constant h = 6.626 10^-34 Js).

The partition function provides the link with key thermodynamic functions. In the case of the canonical ensemble, the Helmholtz free energy A is expressed as:

\[A = -k_{B}T ln Q_{NVT}\]

Using this relationship, numerous properties can be derived. [3] If we have a finite collection of system states, which is representative of the NVT statistical ensemble, average properties can be obtained by simple arithmetic averaging over the n states composing the collection without explicitly computing the partition function. For instance, the average energy is:

(1)\[\langle E \rangle = \dfrac{1}{n} \sum_{i=1}^{n} E_{i}\]

Other statistical ensembles can be defined in the same way as the canonical ensemble as summarized below:

• If an intensive variable is fixed, the associated extensive variable fluctuates.
• If temperature is fixed, energy fluctuates.
• If pressure is fixed, volume fluctuates,
• If chemical potential is fixed, the corresponding number of molecules fluctuates.

The NVT and NPT ensembles are available in MedeA LAMMPS and MedeA GIBBS, but the Grand Canonical Ensemble and the Gibbs Ensemble are available from MedeA GIBBS only)

Statistical ensemble	Imposed variables	Associated thermodynamic potential	Applications
Canonical ensemble	N,V,T	A = E - TS (Helmholtz free) energy)	Phase properties (P,H,Cv,μv…)
Grand Canonical Ensemble	μi,V,T	PV (i.e. \(E-TS-\sum\mu_{i}N_{i}\) )	Adsorption in microporous solids
Isothermal - isobaric ensemble	N,P,T	G = H - TS (Gibbs free energy)	Phase properties \(H, C_{p}, \rho, \mu_{v}, ...\)
Gibbs Ensemble at imposed global volume (m phases)	\(N = N_{1} + ... + N_{m}\) \(V = V_{1} + ... + V_{m}\) T	A = E - TS (Helmholtz free energy)	Phase equilibrium of pure components and mixtures
Gibbs Ensemble at imposed pressure (m phases)	\(N = N_{1} + ... + N_{m}\), P, T	G = H - TS (Gibbs free energy)	Phase equilibrium of mixtures

3.3.2. NPT ensemble

If we want to simulate a system at imposed pressure and temperature, we use the isothermal-isobaric ensemble or NPT ensemble, where volume and energy are fluctuating variables. The probability of any given state-with total energy \(E_{j}\) and volume \(V_{j}\) - is proportional to:

\[exp \left( -\beta E_{j} -\beta PV \right)\]

This ensemble may be used to compute the average energy using Eq. (1) and the average volume using:

\[\langle V \rangle = \dfrac{1}{n} \sum_{j=1}^{n} V_{j}\]

3.3.2.1. Grand Canonical or \({\mu}\)VT Ensemble

If we want to simulate adsorption isotherms, either for pure compounds or multicomponent systems, we must be able to impose temperature and partial pressures to simulate a given point on the isotherm. For this purpose, the most convenient statistical ensemble is the grand canonical ensemble, also termed the \({\mu}\)VT ensemble. At low pressure, the chemical potential \(\mu_{i}\) is linked with the partial pressures \(P_{i} = Py_{i}\) where P is total gas pressure \(y_{i}\) is the molar fraction of component i. A more general expression involves the fugacities \(f_{i}\) in the fluid phase in equilibrium with the adsorbent, which are equivalent to partial pressure in the low-pressure limit:

\[\mu_{i} = \mu_{i_{0}} + RTln\dfrac{f_{i}}{P_{0}} \approx \mu_{i_{0}} + RT ln \dfrac{y_{i}P}{P_{0}}\]

where \(P_{0}\) is a reference pressure. The chemical potential of each species appears explicitly in the probability of each system state in the grand canonical ensemble, which is proportional to:

\[exp \left( -\beta E_{j} + \beta N_{1} \mu_{1} + \beta N_{2}\mu_{2} + ... \right)\]

A grand canonical simulation provides the average number of molecules \(\langle N_{i} \rangle\) in the system for every species i. In the case of microporous adsorbents, this is the average number of adsorbed molecules. The heat of adsorption may be also computed.

3.3.2.2. Microcanonical or NVE Ensemble

The microcanonical or NVE ensemble, in which the number of molecules, global volume, and internal energy are fixed, is obtained by solving the equations of motion in molecular dynamics without thermostat or barostat. In this ensemble, every system state has the same probability. Temperature fluctuates, and the average temperature can be obtained by considering the average kinetic energy per degree of freedom of the system:

\[T = \dfrac{2 \langle K \rangle}{k_{B}N_{f}}\]

where K is the kinetic energy of the system, expressed as a function of the masses \(m_{i}\) and the velocities \(v_{i}\) of its N constituent particles:

\[K = \sum_{i=1}^{N} \dfrac{m_{i}v_{i}^{2}}{2}\]

and \(N_{f} = 3N - N_{c}\) is the total number of degrees of freedom, \(N_{c}\) being the total number of independent constraints such as imposed bond lengths or imposed bond angles.

3.3.2.3. Gibbs Ensemble

The explicit simulation of liquid-vapor interfaces, using either the NVT or NPT statistical ensemble, requires large systems. It is mostly used to investigate the properties of the interface itself. A more efficient way of computing phase equilibria by molecular simulation uses the Gibbs Ensemble [4], where a separate simulation box is used for each phase without any explicit interface. In this ensemble, both temperature and the total number of molecules are fixed, and we can impose either global volume (i.e. the sum of phase volumes) or pressure. As a result of transfer moves from one phase to the other, the number of molecules displays opposite fluctuations in each phase, and chemical potentials are equal in both phases. The sketch shows a two-phase Gibbs ensemble for a system comprising three molecular types (cyclohexane, n-butane, and methane). The arrows illustrate the transfer of molecules from one box to the other, which is specific of this statistical ensemble.

When pressure is imposed in the Gibbs ensemble, phase volumes fluctuate independently. Because of the phase rule that limits the number of imposed intensive parameters, this option can be used only when more than one component is considered.

When global volume is imposed in the Gibbs ensemble, phase volumes fluctuate in opposite directions. This option is applicable either to single or multicomponent systems. In the case of a pure component, it allows the determination of its vapor-liquid coexistence properties.

The Gibbs ensemble is a particular case of the NVT or NPT statistical ensemble in which interface energy is not accounted for.

Phase volumes and energies may be derived by applying the usual averaging procedure defined by Eq. (1). Pressure may be computed for each phase from the virial equation (see next section). This may be used to check that mechanical equilibrium is reached because pressure must be equal in both phases (to within the statistical uncertainties).

3.3.3. Determination of Average Properties

Average pressure is obtained by the Virial equation:

\[\langle P \rangle = \dfrac{Nk_{B}T}{V} + \langle W_{j} \rangle\]

where \(W_{j}\) is a term called the Virial, which accounts for intermolecular interactions through the forces acting on molecule k in the j^th configuration of the ensemble:

\[\langle W_{j} \rangle = \dfrac{1}{3V} \sum r_{k} \cdot F_{k}\]

where \(r_{k}\) is the position of the molecular center of mass in configuration j. If we recall that \(k_{B} = R/N_{A}\), where \(N_{A}\) is Avogadro’s number, the expression reduces to the perfect gas law (PV = RT for one mole) at very low densities when intermolecular forces tend to zero.

The average density (in molecules per unit volume) is derived using

\[\rho = \dfrac{N}{\langle V \rangle}\]

This is different from the arithmetic average of densities, which has no physical sense.

In statistical ensembles where temperature is imposed, deriving the average energy of the system through equation (1) is straightforward. The reference state for energy corresponds to zero kinetic energy and zero potential energy, i.e. the conditions of a dilute gas at 0 K. Once energy is known, enthalpy can be derived through \(\langle H \rangle = \langle E \rangle + \langle P \rangle V\) in the canonical ensemble or \(\langle H \rangle = \langle E \rangle + P \langle V \rangle\) in the NPT ensemble.

The chemical potential is defined as a derivative of the Helmholtz energy:

\[\mu_{i} \left( \dfrac{\partial A}{\partial N_{i}} \right)_{V,T,N_{j, j \ne i}}\]

However, it cannot be evaluated simply as an ensemble average.

For extensive properties, comparison of simulation results with macroscopic measurements is more convenient when properties are expressed on a molar basis. For instance, the molar volume v may be computed from an NPT simulation using:

\[v = \dfrac{N_{a}}{N} \langle V \rangle\]

The same conversion may be applied to other extensive variables like energy and enthalpy. Intensive variables such as pressure or chemical potential do not need any correction factor to be compared with measurements, provided a consistent unit system is used.

3.3.4. Determination of Derivative Properties from Fluctuations

The amplitude of the fluctuations around the average makes it possible to determine partial derivatives of the property considered. For instance, the standard deviation of energy in the canonical ensemble is linked with the partial derivative of energy versus temperature, i.e. the heat capacity \(C_{V}\) of the system:

\[C_{V} = \left( \dfrac{\partial E}{\partial T} \right) _{V} = \dfrac{1}{k_{B}T^{2}} \left( \langle E^{2} \rangle - \langle E \rangle ^{2} \right)\]

In this expression, the right-hand side represents energy fluctuations around the mean value, which can be computed by molecular simulation. Similarly, fluctuations in volume in the NPT ensemble may be used to obtain the compressibility coefficient:

\[\beta_{T} = - \dfrac{1}{V} \left( \dfrac{\partial V}{\partial P} \right) _{T} = \dfrac{1}{k_{B}T^{2}} \left( \langle V^{2} \rangle - \langle V \rangle ^{2} \right)\]

Comparable fluctuation formulae are available for other derivative properties.

3.3.5. Possible Ways of Simulating Ensembles

The simplest way to build a statistical ensemble is to solve the equations of motion for every atom in every molecule, accounting for intramolecular and intermolecular interactions inside the system: this way is known as molecular dynamics. An alternate way, known as Monte Carlo simulation, consists of using statistical methods to generate the collection of configurations with the desired probability distribution without solving the equations. Both methods yield identical ensemble averages (this is known as the ergodicity theorem).

The Newton equations of motion, which have to be integrated in molecular dynamics, are expressed in the form:

\[F_{i} = m_{i} \dfrac{d^{2}r_{i}}{dt^{2}}\]

where r_i is the position of an atom, i.e. a three dimensional vector (\(x_{i}\), \(y_{i}\), \(z_{i}\))

\(F_{i}\) is the force acting on the atom from the other atoms of the system, which :sub:`` can be computed by deriving the potential energy U versus coordinates:

\[F_{i} = \nabla U (r_{i})\]

The solution of the equations of motion is such that total energy (i.e. the sum of kinetic and potential energy) is constant. Therefore, molecular dynamics simulates the microcanonical or NVE ensemble, unless specific steps are taken to force energy or volume to fluctuate.

In order to perform a molecular dynamics calculation at imposed temperature, the kinetic energy is changed by altering the center of mass velocities at regular intervals: this is the purpose of thermostat techniques, which allow simulation in the NVT ensemble [5].

Molecular dynamics at imposed pressure is possible through barostat or barostat methods that allows for volume variations [6] to obtain the correct probability distribution equation.

It is difficult to simulate open ensembles (Grand Canonical and Gibbs ensemble) with molecular dynamics because molecule insertions or deletions cause discontinuities in the integration of Newton’s equations. Monte Carlo is more often used for such ensembles.

Molecular dynamics is most useful for simulating the dynamic behavior of a system, particularly its transport properties. By transport properties, we mean here the coefficients governing the rate of mass transfer (diffusion coefficients), heat transfer (thermal conductivity) or momentum transfer (viscosity). For instance, the self-diffusion coefficient of a component in a space of dimension d can be expressed from the mean squared displacement, according to Einstein’s expression:

\[D = \dfrac{\langle \left( r(t) - r(t_{0}) \right)^{2} \rangle}{2d \left( t - t_{0} \right)}\]

Shear viscosity \({\eta}\) and thermal conductivity \({\lambda}\) may be also determined using molecular dynamics through Green Kubo expressions (see Ungerer, Nietro-Draghi, Rousseau, Ahunbay & Lachet [7] and references therein):

\[\eta = \dfrac{V}{k_{B}T} \int _{0} ^{\inf} \langle \sigma_{xy}(t) \sigma_{xy}(0) \rangle\]

\[\lambda = \dfrac{V}{3k_{B}T} \int _{0} ^{\inf} \langle J_{q}(t) J_{q}(0) \rangle dt\]

where \(\sigma_{xy}(t)\) and \(J_{q}(t)\) are pressure tensor components and heat flux. It is also possible to use non-equilibrium molecular dynamics, in which perturbations or specific boundary conditions are imposed, to get transport coefficients [8].

[1]	Philippe Ungerer, B Tavitian, and A. Boutin, Applications of Molecular Simulation in the Oil and Gas Industry: Monte Carlo Methods, (Editions Technip, 2005).

[2]

D.A. McQuarrie, Statistical Mechanics, (Harper & Row, 1975); M P Allen and D J Tildesley, Computer Simulation of Liquids, Book, (Clarendon Press, 1987); Erwin Schrödinger, Statistical Thermodynamics, (Dover Publications, 1989); Daan Frenkel and Berend Smit, Understanding Molecular Simulation: From Algorithms to Applications, 1st ed., (Academic Press, 1996).

[3]	D.A. McQuarrie, Statistical Mechanics, (Harper & Row, 1975).

[4]	Athanassios Panagiotopoulos, “Direct Determination of Phase Coexistence Properties of Fluids by Monte Carlo Simulation in a New Ensemble,” Molecular Physics 61, no. 4 (July 1, 1987): 813-826.

[5]	Shuichi Nos.. é, “A Unified Formulation of the Constant Temperature Molecular Dynamics Methods,” Journal of Chemical Physics 81, no. 1 (1984): 511. William Hoover, “Canonical Dynamics: Equilibrium Phase-Space Distributions,” Physical Review A 31, no. 3 (March 1985): 1695-1697.

[6]	D Evans and G Morriss, “Isothermal-Isobaric Molecular Dynamics,” Chemical Physics 77, no. 1 (May 15, 1983): 63-66.

[7]	Philippe Ungerer, Carlos Nieto-Draghi, Bernard Rousseau, Göktug Ahunbay, and Véronique Lachet, “Molecular Simulation of the Thermophysical Properties of Fluids: From Understanding Toward Quantitative Predictions,” Journal of Molecular Liquids 134, no. 1 (May 15, 2007): 71-89.

[8]

B Y Wang and P T Cummings, “Non-Equilibrium Molecular Dynamics Calculation of the Shear Viscosity of Carbon Dioxide/Ethane Mixtures,” Molecular Simulation 10, no. 1 (April 1, 1993): 1-11; RL Rowley, Statistical Mechanics for Thermophysical Property Prediction, Prentice Hall, 1994; F Müller-Plathe and D Reith, “Cause and Effect Reversed in Non-Equilibrium Molecular Dynamics: an Easy Route to Transport Coefficients,” Computational and Theoretical Polymer Science 9, no. 3 (1999): 203-209.

download:pdf