3.3. Statistical Thermodynamics

---

download:pdf

---

This chapter’s purpose, largely inspired by a textbook by Ungerer, Tavitian & Boutin [1], is to help MedeA GIBBS and MedeA LAMMPS users 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 complex treatment is generally unnecessary for thermodynamic properties, as the related molecular interactions do not involve significant electronic structure modifications. In compensation, molecular simulation can consider larger systems which 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 the section Statistical ensembles and partition functions. As the distribution of energy among molecules plays a central role in statistical mechanics, we will then review how 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 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 for computing 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, of dimension 6N for N particles, is called the phase space.

There are several different types of statistical ensemble, depending on the type of physical 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 called the canonical 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:

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

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

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

Where

(3)\[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 corresponding to the same state. Consistently with quantum mechanics, the finite summation may be transformed in a continuous integral:

(4)\[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:

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

Using this relationship, we can derive numerous properties. [3] If we have a finite collection of system states, 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:

(6)\[\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 the temperature is fixed, energy fluctuates.
• If the pressure is fixed, volume fluctuates,
• If the chemical potential is fixed, the corresponding number of molecules fluctuates.

The NVT and NPT ensembles are available in both MedeA LAMMPS and MedeA GIBBS, but the Grand Canonical Ensemble, the Osmotic Ensemble, the Semi-grand Canonical Ensemble, and the Gibbs Ensemble are available in 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:

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

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

(8)\[\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 pressures in the low-pressure limit:

(9)\[\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:

(10)\[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 also be 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 a thermostat or a barostat. In this ensemble, every system state has the same probability. The temperature fluctuates, and the average value can be obtained by considering the average kinetic energy per degree of freedom of the system:

(11)\[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:

(12)\[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, and it is mainly 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 to this statistical ensemble.

When pressure is imposed in the Gibbs ensemble, phase volumes fluctuate independently. Because the phase rule limits the number of imposed intensiveparameters, 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. (6). Pressure may be computed for each phase from the virial equation (see next section). The computation of pressure may be used to check that mechanical equilibrium is reached because, at equlibrium, pressure must be equal in both phases (within the statistical uncertainties).

3.3.3. Determination of Average Properties

The Virial equation provides the average pressure:

(13)\[\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:

(14)\[\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

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

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

In statistical ensembles where the temperature is imposed, deriving the average energy of the system through equation (6) 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:

(16)\[\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:

(17)\[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:

(18)\[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:

(19)\[\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, including the isosteric heat of adsorption.

3.3.4.1. Isosteric Heat of adsorption

The differential molar enthalpy of adsorption (ΔΗ^s), also called the isosteric enthalpy of adsorption (Q_st), is defined as:

(20)\[\Delta H = H_{\text{sorbed}} - H_{\text{gas}}\]

where \(H_{\text{sorbed}}\) and \(H_{\text{gas}}\) are the partial molar enthalpies in the sorbed and gas phases, respectively. The enthalpy is the sum of the internal energy and the product of pressure and volume, PV. For the vapor phase, PV is assumed to be equal to RT, and the adsorbed phase’s molecular volume is neglected. Assuming molar kinetic energy is equal in the gas and the adsorbed state, we can express the heat of adsorption as a function of the total molar potential energy in the gas phase and the adsorbed phase:

(21)\[-\Delta H^{0} = RT - \overline{E}_{\text{sorbed}}^{\text{tot}} + \overline{E}_{\text{gas}}^{\text{tot}}\]

In GCMC simulations, it is equivalent to calculating -\({\Delta H}^{0}\) using the partial derivative of the average total energy with respect to the average number of adsorbed molecules in both gas and adsorbed phases:

(22)\[\Delta H^{0} = RT - \frac{\partial\left\langle E_{\text{sorbed}}^{\text{tot}} \right\rangle}{\partial\left\langle N_{\text{sorbed}} \right\rangle} + \frac{\partial\left\langle E_{\text{gas}}^{\text{tot}} \right\rangle}{\partial\left\langle N_{\text{gas}} \right\rangle}\]

where \(\left\langle N_{\text{sorbed}} \right\rangle\) and \(\left\langle N_{\text{gas}} \right\rangle\) are the average number of molecules in the adsorbed and gas phases, respectively. The heat of adsorption \(Q_{\text{st}}\) can thus be calculated by the fluctuations method:

(23)\[Q_{\text{st}} = RT - \frac{\left\langle E_{\text{sorbed}}^{\text{tot}}N \right\rangle - \left\langle E_{\text{sorbed}}^{\text{tot}} \right\rangle\left\langle N \right\rangle}{\left\langle N_{\text{sorbed}}^{2} \right\rangle - \left\langle N_{\text{sorbed}} \right\rangle^{2}} + \frac{\left\langle E_{\text{gas}}^{\text{tot}}N \right\rangle - \left\langle E_{\text{gas}}^{\text{tot}} \right\rangle\left\langle N \right\rangle}{\left\langle N_{\text{gas}}^{2} \right\rangle - \left\langle N_{\text{gas}} \right\rangle^{2}}\]

The heat of adsorption calculation requires a Monte Carlo simulation of the vapor phase in the Grand Canonical ensemble to determine the last part of the above equation ((23)). These simulations are usually performed so that the number of molecules fluctuates reasonably around five or six, which often implies simulation box dimensions of about 1,000 Å. However, usually the calculation of \(Q_{\text{st}}\) is achieved without this second simulation of the vapor phase.

If the vapor phase is assumed to be ideal, the third part of the equation is equivalent to the molar intramolecular energy of molecules in the vapor phase. Moreover, if internal degrees of freedom are considered not to be affected by adsorption, the molar intramolecular energies of the vapor phase and the adsorbed phase are equal. These approximations lead to the following expression for the heat of adsorption, proposed by Nicholson and Parsonage:

(24)\[Q_{\text{st}} = RT - \frac{\left\langle E_{\text{sorbed}}^{\text{ext}}N \right\rangle - \left\langle E_{\text{sorbed}}^{\text{ext}} \right\rangle\left\langle N \right\rangle}{\left\langle N_{}^{2} \right\rangle - \left\langle N_{} \right\rangle^{2}}\]

where \(\left\langle E_{\text{sorbed}}^{\text{ext}} \right\rangle\) represents intermolecular interactions in the adsorbed phase.

3.3.4.2. Grid Isosteric Heat of Adsorption

Eq. (24), commondly used to determine the isosteric heat of adsorption, is based on severe approximations: ideal gas for the fluid phase and the assumption that intramolecular energies in the vapor phase and the adsorbed are equal. This last hypothesis is only sometimes fulfilled. Nevertheless, at low coverage, when sorbate/sorbate interactions are feeble, the isosteric heat \(Q_{\text{st}}\) can be determined through equation (24) by considering only the adsorbate/zeolite interactions in the intermolecular potential energy (see (25)):

(25)\[Q_{st,0} = RT - \frac{\left\langle E_{\text{sobed}}^{\text{grid}}N \right\rangle - \left\langle E_{\text{sorbed}}^{\text{grid}} \right\rangle\left\langle N \right\rangle}{\left\langle N_{}^{2} \right\rangle - \left\langle N_{} \right\rangle^{2}}\]

where \(E_{\text{sorbed}}^{\text{grid}}\) is the interaction energy of adsorbate and adsorbent (i.e., neglecting adsorbate-adsorbate interactions).

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).

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

(26)\[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:

(27)\[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.

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 allow for volume variations [6] to obtain the correct probability distribution equation.

Simulating open ensembles (Grand Canonical and Gibbs ensemble) with molecular dynamics is difficult because molecule insertions or deletions cause discontinuities in integrating 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:

(28)\[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 also be determined using molecular dynamics through Green Kubo expressions (see Ungerer, Nietro-Draghi, Rousseau, Ahunbay & Lachet [7] and references therein):

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

(30)\[\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