Editing Gibbs-Duhem integration

The so-called '''Gibbs-Duhem integration''' refers  to a number of methods that couple
molecular [[Computer simulation techniques |simulation techniques]]  with [[Thermodynamic relations |thermodynamic equations]] in order to draw
[[Computation of phase equilibria | phase coexistence]] lines. The original method was proposed by David Kofke <ref>[http://dx.doi.org/10.1080/00268979300100881 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)]</ref>
<ref>[http://dx.doi.org/10.1063/1.465023 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)]</ref>.

== Basic Features ==

Consider two thermodynamic phases: <math> a </math> and <math> b  </math>,  at thermodynamic equilibrium at certain conditions. Thermodynamic equilibrium implies:

* Equal [[temperature]] in both phases: <math> T = T_{a} = T_{b} </math>, i.e. thermal equilibrium.
* Equal [[pressure]] in both phases <math> p = p_{a} = p_{b} </math>, i.e. mechanical equilibrium.
* Equal [[chemical potential]]s for the components <math> \mu_i = \mu_{ia} = \mu_{ib} </math>, i.e. ''material'' equilibrium.

In addition, if one is  dealing with a statistical mechanical [[models |model]], having certain parameters that can be represented as <math> \lambda </math>, 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 <math> T , p, \lambda </math> two phases of the systems are at equilibrium, this implies:

: <math> \mu_{a} \left( T, p, \lambda \right) = \mu_{b} \left( T, p, \lambda \right) </math>

Given the thermal equilibrium we can also write:

: <math> \beta \mu_{a} \left( \beta, \beta p, \lambda \right) = \beta \mu_{b} \left( \beta, \beta p, \lambda \right) </math>

where
* <math> \beta := 1/k_B T </math>, where <math> k_B </math> is the [[Boltzmann constant]]
When a differential change of the conditions is performed one will have, for any phase:

: <math> d \left( \beta\mu \right) = \left[ \frac{ \partial (\beta \mu) }{\partial \beta} \right]_{\beta p,\lambda} d \beta +
\left[ \frac{ \partial (\beta \mu) }{\partial (\beta p)} \right]_{\beta,\lambda} d (\beta p) + 
\left[ \frac{ \partial (\beta \mu) }{\partial \lambda} \right]_{\beta,\beta p} d \lambda.
</math> 

Taking into account that <math> \mu </math> is the [[Gibbs energy function]] per particle
: <math> d \left( \beta\mu \right) =  \frac{E}{N}  d \beta +  \frac{ V }{N } d (\beta p)  + 
\left[ \frac{ \partial (\beta \mu) }{\partial \lambda} \right]_{\beta,\beta p} d \lambda.
</math> 

where:
* <math> \left. E \right. </math> is the [[internal energy]] (sometimes written as <math>U</math>).

* <math> \left. V \right. </math> is the volume

* <math> \left. N \right. </math> is the number of particles

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

: <math> \bar{L} \equiv \frac{1}{N} \left[ \frac {\partial (\beta \mu )}{\partial \lambda }\right]_{\beta,\beta p} </math>

Again, let us suppose that we have a phase coexistence at a point given by <math>\left[ \beta_0, (\beta p)_0, \lambda_0 \right]</math> and that
we want to modify slightly the conditions. In order to keep the system at the coexistence conditions:

: <math> d \left[ \beta \mu_{a} - \beta \mu_b \right] = 0 </math>

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

:<math>  \left( \Delta  \bar{E} \right) d  \beta + \left( \Delta \bar{V} \right) d (\beta p) + \left(\Delta  \bar{L} \right) d \lambda = 0 </math>

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

:<math> d(\beta p) = - \frac{ \Delta \bar{L} }{\Delta \bar{V} } d \lambda. </math> (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 [[Isothermal-isobaric ensemble |NpT ensemble]]) runs to estimate the values of <math> \bar{L}, \bar{V} </math> for both
phases at given values of <math> [\beta, \beta p,  \lambda ] </math>.

* 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>[http://dx.doi.org/10.1063/1.2137705      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)]</ref> and [[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 ==
<references/>
'''Related reading'''
*[http://dx.doi.org/10.1063/1.2137706      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)]
*[http://dx.doi.org/10.1063/1.3486090  Gerassimos Orkoulas "Communication: Tracing phase boundaries via molecular simulation: An alternative to the Gibbs–Duhem integration method", Journal of Chemical Physics '''133''' 111104 (2010)]

[[category: computer simulation techniques]]