Gibbs-Duhem integration

The so-called Gibbs-Duhem integration refers to a number of methods that couple molecular simulation techniques with thermodynamic equations in order to draw phase coexistence lines. The original method was proposed by David Kofke (Refs. 1 and 2).

Basic Features

Consider two thermodynamic phases: $a$ and $b$, at thermodynamic equilibrium at certain conditions. Thermodynamic equilibrium implies:

• Equal temperature in both phases: $T=T_a=T_b$, i.e. thermal equilibrium.
• Equal pressure in both phases $p=p_a=p_b$, i.e. mechanical equilibrium.
• Equal chemical potentials for the components ${\mu }_i={\mu }_{ia}={\mu }_{ib}$, i.e. material equilibrium.

In addition, if one is dealing with a statistical mechanical model, having certain parameters that can be represented as $\lambda$, then the model should be the same in both phases.

Example: phase equilibria of one-component system

Notice: The derivation that follows is just a particular route to perform the integration

• Consider that at given conditions of $T,p,\lambda$ two phases of the systems are at equilibrium, this implies:

${\mu }_a(T,p,\lambda )={\mu }_b(T,p,\lambda )$

Given the thermal equilibrium we can also write:

$\beta {\mu }_a(\beta ,\beta p,\lambda )=\beta {\mu }_b(\beta ,\beta p,\lambda )$

where

• $\beta :=1/k_{B}T$, where $k_B$ is the Boltzmann constant

When a differential change of the conditions is performed one will, have for any phase:

$d(\beta \mu )={[\frac{\partial (\beta \mu )}{\partial \beta }]}_{\beta p,\lambda }d\beta +{[\frac{\partial (\beta \mu )}{\partial (\beta p)}]}_{\beta ,\lambda }d(\beta p)+{[\frac{\partial (\beta \mu )}{\partial \lambda }]}_{\beta ,\beta p}d\lambda .$

Taking into account that $\mu$ is the Gibbs energy function per particle

$d(\beta \mu )=\frac{E}{N}d\beta +\frac{V}{N}d(\beta p)+{[\frac{\partial (\beta \mu )}{\partial \lambda }]}_{\beta ,\beta p}d\lambda .$

where:

• $E$ is the internal energy (sometimes written as $U$).

• $V$ is the volume

• $N$ is the number of particles

$E,V$ are the mean values of the energy and volume for a system of $N$ particles in the isothermal-isobaric ensemble

Let us use a bar to design quantities divided by the number of particles: e.g. $\overline{E}=E/N;\overline{V}=V/N$; and taking into account the definition:

$\overline{L}\equiv {[\frac{\partial (\beta \mu )}{\partial \lambda }]}_{\beta ,\beta p}$

Again, let us suppose that we have a phase coexistence at a point given by $[{\beta }_0,(\beta p{)}_0,{\lambda }_0]$ and that we want to modify slightly the conditions. In order to keep the system at the coexistence conditions:

$d[\beta {\mu }_a-\beta {\mu }_b]=0$

Therefore, to keep the system on the coexistence conditions, the changes in the variables $\beta ,(\beta p),\lambda$ are constrained to fulfill:

$(\mathrm{\Delta }\overline{E})d\beta +(\mathrm{\Delta }\overline{V})d(\beta p)+(\mathrm{\Delta }\overline{L})d\lambda =0$

where for any property $X$ we can define: $\mathrm{\Delta }X\equiv X_a-X_b$ (i.e. the difference between the values of the property in the phases). Taking a path with, for instance constant $\beta$, the coexistence line will follow the trajectory produced by the solution of the differential equation:

$d(\beta p)=-\frac{\mathrm{\Delta }\overline{L}}{\mathrm{\Delta }\overline{V}}d\lambda .$ (Eq. 1)

The Gibbs-Duhem integration technique, for this example, will be a numerical procedure covering the following tasks:

• Computer simulation (for instance using Metropolis Monte Carlo in the NpT ensemble) runs to estimate the values of $\overline{L},\overline{V}$ for both

phases at given values of $[\beta ,\beta p,\lambda ]$.

• A procedure to solve numerically the differential equation (Eq.1)

Peculiarities of the method (Warnings)

• A good initial point must be known to start the procedure (See Ref. 3 and the entry: computation of phase equilibria).

• The integrand of the differential equation is computed with some numerical uncertainty

• Care must be taken to reduce (and estimate) possible departures from the correct coexistence lines

References

1. David A. Kofke, "Gibbs-Duhem integration: a new method for direct evaluation of phase coexistence by molecular simulation", Molecular Physics 78 pp 1331 - 1336 (1993)
2. David A. Kofke, "Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line", Journal of Chemical Physics 98 pp. 4149-4162 (1993)
3. A. van 't Hof, S. W. de Leeuw, and C. J. Peters "Computing the starting state for Gibbs-Duhem integration", Journal of Chemical Physics 124 054905 (2006)
4. A. van 't Hof, C. J. Peters, and S. W. de Leeuw "An advanced Gibbs-Duhem integration method: Theory and applications", Journal of Chemical Physics 124 054906 (2006)