Editing Computation of phase equilibria
Jump to navigation
Jump to search
The edit can be undone. Please check the comparison below to verify that this is what you want to do, and then publish the changes below to finish undoing the edit.
Latest revision | Your text | ||
Line 1: | Line 1: | ||
Thermodynamic equilibrium implies, for two phases <math> \alpha </math> and <math> \beta </math>: | Thermodynamic equilibrium implies, for two phases <math> \alpha </math> and <math> \beta </math>: | ||
* | |||
* | * Equal [[temperature]]: <math> T_{\alpha} = T_{\beta} </math> | ||
* | |||
* Equal [[pressure]]: <math> p_{\alpha} = p_{\beta} </math> | |||
* Equal [[chemical potential]]: <math> \mu_{\alpha} = \mu_{\beta} </math> | |||
The computation of phase equilibria using computer simulation can follow a number of different strategies. Here we will focus mainly | |||
on first order transitions in fluid phases, usually liquid-vapor equilibria. | |||
== Independent simulations for each phase at fixed temperature in the [[canonical ensemble]] == | == Independent simulations for each phase at fixed temperature in the [[canonical ensemble]] == | ||
Simulations can be carried out using either the [[Monte Carlo]] or the [[molecular dynamics]] technique. | Simulations can be carried out using either the [[Monte Carlo]] or the [[molecular dynamics]] technique. | ||
Assuming that one has some knowledge on the | Assuming that one has some knowledge on the phase diagram of the system, one can try the following recipe: | ||
# Fix a temperature and a number of particles | # Fix a temperature and a number of particles | ||
#Perform a limited number of simulations in the low density region (where the gas phase density is expected to be) | #Perform a limited number of simulations in the low density region (where the gas phase density is expected to be) | ||
Line 20: | Line 25: | ||
because is is not unusual have large uncertainties in the results for the properties. | because is is not unusual have large uncertainties in the results for the properties. | ||
The basic idea is to use [[thermodynamic consistency]] requirements to improve the analysis. | The basic idea is to use [[thermodynamic consistency]] requirements to improve the analysis. | ||
== Methodology in the [[Isothermal-isobaric ensemble|NpT]] ensemble == | == Methodology in the [[Isothermal-isobaric ensemble|NpT]] ensemble == | ||
=== Low temperature: <math> \left. T << T_c \right. </math> === | === Low temperature: <math> \left. T << T_c \right. </math> === | ||
For temperatures well below the | |||
For temperatures well below the critical point, provided that the calculation of the chemical potential | |||
of the liquid phase using [[Widom test-particle method]] gives precise results, the following strategy can be used to obtain a 'quick' result: | of the liquid phase using [[Widom test-particle method]] gives precise results, the following strategy can be used to obtain a 'quick' result: | ||
#Perform an <math> NpT </math> simulation of the liquid phase at zero pressure, i.e. <math> p \simeq 0 </math> | #Perform an <math> NpT </math> simulation of the liquid phase at zero pressure, i.e. <math> p \simeq 0 </math> | ||
#Arrive at an initial estimate, <math> \mu^{(1)} </math> for the coexistence value of the chemical potential by computing, in the liquid phase: <math> \left. \mu^{(1)} = \mu_l (N,T,p=0) \right. </math> | #Arrive at an initial estimate, <math> \mu^{(1)} </math> for the coexistence value of the chemical potential by computing, in the liquid phase: <math> \left. \mu^{(1)} = \mu_l (N,T,p=0) \right. </math> | ||
#Make a first estimate of the coexistence pressure, <math> p^{(1)} </math>, by computing, either via simulation or via the [[Virial coefficients of model systems |virial coefficients]] of the gas phase, the pressure at which the gas phase | #Make a first estimate of the coexistence pressure, <math> p^{(1)} </math>, by computing, either via simulation or via the [[Virial coefficients of model systems |virial coefficients]] of the gas phase, the pressure at which the gas phase fulfills: <math> \left. \mu_g(N,T, p^{(1)} ) = \mu^{(1)} \right. </math> | ||
#Refine the results, if required, by performing a simulation of the liquid phase at <math> \left. p^{(1)} \right. </math>, or use estimates of <math> \left( \partial V / \partial p \right)_{N,T,p=0} </math> (from the initial simulation) and the gas equation of state data to correct the initial estimates of pressure and chemical potential at coexistence. | #Refine the results, if required, by performing a simulation of the liquid phase at <math> \left. p^{(1)} \right. </math>, or use estimates of <math> \left( \partial V / \partial p \right)_{N,T,p=0} </math> (from the initial simulation) and the gas equation of state data to correct the initial estimates of pressure and chemical potential at coexistence. | ||
Note that this method works only if the liquid phase remains metastable at zero pressure. | Note that this method works only if the liquid phase remains metastable at zero pressure. | ||
=== Weak first order transitions === | === Weak first order transitions === | ||
There are situations where other strategies based on [[Isothermal-isobaric ensemble|NpT]] ensemble simulations can be used. | There are situations where other strategies based on [[Isothermal-isobaric ensemble|NpT]] ensemble simulations can be used. | ||
Similar approaches have also been applied in the [[Monte Carlo in the grand-canonical ensemble|Grand Canonical]] ensemble. | Similar approaches have also been applied in the [[Monte Carlo in the grand-canonical ensemble|Grand Canonical]] ensemble. | ||
Line 35: | Line 45: | ||
to cross from one phase to the other when it is being simulated. Somehow, the situation is now the opposite to that | to cross from one phase to the other when it is being simulated. Somehow, the situation is now the opposite to that | ||
described in Section 2.1. | described in Section 2.1. | ||
Taking into the account the classical partition function in the [[isothermal-isobaric ensemble|NpT]] ensemble it can be written: | Taking into the account the classical partition function in the [[isothermal-isobaric ensemble|NpT]] ensemble it can be written: | ||
Line 42: | Line 53: | ||
conditions <math> \left. N,P,T \right. </math>. Let <math> p_{eq} </math> the pressure at which the phase transition occurs. In such a | conditions <math> \left. N,P,T \right. </math>. Let <math> p_{eq} </math> the pressure at which the phase transition occurs. In such a | ||
case the following scenario is expected for <math> \left. P(V|N,p,T) \right. </math>: | case the following scenario is expected for <math> \left. P(V|N,p,T) \right. </math>: | ||
*<math> \left. P(V|N,p_{eq},T) \right. </math> has two maxima, corresponding to the liquid and vapor pure phases, with <math> \left. P(V_v|N,p_{eq},T) = P(V_l|N,p_{eq},T) = P_{v/l} \right. </math> | *<math> \left. P(V|N,p_{eq},T) \right. </math> has two maxima, corresponding to the liquid and vapor pure phases, with | ||
:: <math> \left. P(V_v|N,p_{eq},T) = P(V_l|N,p_{eq},T) = P_{v/l} \right. </math> | |||
*The probability of a given intermediate volume at <math> \left. p_{eq} \right. </math> can be estimated (from macroscopic arguments) as: | *The probability of a given intermediate volume at <math> \left. p_{eq} \right. </math> can be estimated (from macroscopic arguments) as: | ||
:<math> \left. P(V|N,p_{eq},T) \simeq P_{v/l} \times \exp \left[ - \frac{ \gamma(T) \mathcal A }{k_B T } \right] \right. </math>, | :: <math> \left. P(V|N,p_{eq},T) \simeq P_{v/l} \times \exp \left[ - \frac{ \gamma(T) \mathcal A }{k_B T } \right] \right. </math>, | ||
where <math> \left. \gamma(T) \right. </math> is the [[surface tension]] of the vapor-liquid interface, | where <math> \left. \gamma(T) \right. </math> is the [[surface tension]] of the vapor-liquid interface, | ||
and <math> \left. {\mathcal A} \right. </math> is the | and <math> \left. {\mathcal A} \right. </math> is the | ||
surface area, which depends on the thermodynamic | surface area, which depends on the thermodynamic | ||
variables <math> \left. (N,V,T) \right. </math> and the | variables <math> \left. (N,V,T) \right. </math> and the geommetry of the simulation box. | ||
For small values of the surface tension, small system sizes and good simulation algorithms it could be possible for | For small values of the surface tension, small system sizes and good simulation algorithms it could be possible for | ||
pressures close to <math> p_{eq} </math> to sample in a simulation the whole region of densities between | pressures close to <math> p_{eq} </math> to sample in a simulation the whole region of densities between | ||
vapor and liquid densities. If such is the case the phase equilbria conditions can be computed | |||
by | by reweighting techniques applied on the volume histogramas of the simulation. | ||
==== Simple | |||
==== Simple reweighting of the volume probability distribution ==== | |||
Suppose that a precise <math> \left. NpT \right. </math> simulation has been carried out at pressure <math> \left. p_0 \right. </math>. From that simulation a volume probability distribution, <math> \left. P_0(V|N,p_0,T) \right. </math> | Suppose that a precise <math> \left. NpT \right. </math> simulation has been carried out at pressure <math> \left. p_0 \right. </math>. From that simulation a volume probability distribution, <math> \left. P_0(V|N,p_0,T) \right. </math> | ||
has been computed. It is possible to use this function to estimate the distributions of volume for pressure values | has been computed. It is possible to use this function to estimate the distributions of volume for pressure values | ||
close to <math> \left. p_0 \right. </math>; | close to <math> \left. p_0 \right. </math>; | ||
: <math> \left. P(V|N,p,T) \propto \exp \left[ - \frac{ (p-p_0) V }{k_B T } \right] P_0(V|N,p_0,T) \right. </math> | :: <math> \left. P(V|N,p,T) \propto \exp \left[ - \frac{ (p-p_0) V }{k_B T } \right] P_0(V|N,p_0,T) \right. </math> | ||
Using this | Using this reweighting procedure we can estimate the value of <math> \left. p_{eq} \right. </math>; that is: the value | ||
of pressure for which the volume distribution presents two peaks of equal height. | of pressure for which the volume distribution presents two peaks of equal height. | ||
The analysis of the form of the distributions at the equilibrium conditions for different system sizes (i.e. | The analysis of the form of the distributions at the equilibrium conditions for different system sizes (i.e. | ||
in the current case for different values of <math> \left. N \right. </math>) can be useful to | in the current case for different values of <math> \left. N \right. </math>) can be useful to distinghish between | ||
continuous and first order [[phase transitions]]. | continuous and first order [[phase transitions]]. | ||
[Appropriate references should be quoted here]]--Noe 15:25, 2 October 2007 (CEST) | |||
[This subsection is provisional, I have tho check a number of isssues]--Noe 10:49, 2 October 2007 (CEST) | |||
[Working here] --Noe 10:30, 2 October 2007 (CEST) | |||
=== See also === | === See also === | ||
[[Surface tension]] | [[Surface tension]] | ||
== Van der Waals loops in the [[canonical ensemble|canonical ]] ensemble == | == Van der Waals loops in the [[canonical ensemble|canonical ]] ensemble == | ||
It is possible to compute the liquid- | |||
by performing a number of simulations sampling appropriately the | |||
It is possible to compute the liquid-vapor equilibrium without explicit calculations of the chemical potential (or the pressure) | |||
by performing a number of simulations sampling appropriately the vapor, liquid, and intermediate regions. | |||
As an example, consider a simple fluid at a given subcritical temperature (<math> \left. T < T_c \right. </math>). We can perform a number of simulations for a given number of particles, <math> \left. N \right. </math> and different densities: | As an example, consider a simple fluid at a given subcritical temperature (<math> \left. T < T_c \right. </math>). We can perform a number of simulations for a given number of particles, <math> \left. N \right. </math> and different densities: | ||
Line 78: | Line 105: | ||
In these simulations, we can compute the pressure (or the chemical potential) and fit the result to an appropriate equation. | In these simulations, we can compute the pressure (or the chemical potential) and fit the result to an appropriate equation. | ||
With such an ''equation of state'' the phase | With such an ''equation of state'' the phase equilbria can be estimated. | ||
If two phase equilibria exists, a ''loop'' in the representation of <math> \left. p = p (\rho) \right. </math> (or <math> \left. \mu = \mu (\rho) \right. </math>) | If two phase equilibria exists, a ''loop'' in the representation of <math> \left. p = p (\rho) \right. </math> (or <math> \left. \mu = \mu (\rho) \right. </math>) | ||
should appear. | should appear. | ||
* Computing <math> \left. \mu(\rho) \right. </math> from the equation of state given as <math> \left. p(\rho) \right. </math>: | * Computing <math> \left. \mu(\rho) \right. </math> from the equation of state given as <math> \left. p(\rho) \right. </math>: | ||
For fixed temperature and number of particles: | For fixed temperature and number of particles: | ||
Line 96: | Line 125: | ||
:<math> \left. A = - p V + \mu N ; \; \; \; a = - p/ \rho + \mu \right. </math> | :<math> \left. A = - p V + \mu N ; \; \; \; a = - p/ \rho + \mu \right. </math> | ||
therefore: | |||
therefore: [I have to check the equations]--Noe 13:03, 26 September 2007 (CEST) | |||
: <math> \left. \mu(\rho) - p(\rho)/\rho = \mu(\rho_0) - p(\rho_0)/\rho_0 + \int_{\rho_0}^{\rho} \frac{ p(\rho') }{(\rho')^2} d \rho' \right. </math> | : <math> \left. \mu(\rho) - p(\rho)/\rho = \mu(\rho_0) - p(\rho_0)/\rho_0 + \int_{\rho_0}^{\rho} \frac{ p(\rho') }{(\rho')^2} d \rho' \right. </math> | ||
Line 103: | Line 134: | ||
A similar procedure can be built up to compute <math> \left. p(\rho) \right. </math> | A similar procedure can be built up to compute <math> \left. p(\rho) \right. </math> | ||
from <math> \left. \mu(\rho) \right. </math>. | from <math> \left. \mu(\rho) \right. </math>. | ||
Once <math> \left. p(\rho) \right. </math> and <math> \left. \mu(\rho) \right. </math> are known it is | |||
Once <math> \left. p(\rho) \right. </math> and <math> \left. \mu(\rho) \right. </math> are known it is straighforward to compute the coexistence point. | |||
==== Practical details ==== | ==== Practical details ==== | ||
Some precautions should be taken if this procedure is used: | Some precautions should be taken if this procedure is used: | ||
* The precision of the simulation results in the two phase region will be poor (so, large simulations are required to have a good estimation of the equation of state) | * The precision of the simulation results in the two phase region will be poor (so, large simulations are required to have a good estimation of the equation of state) | ||
* The simulation results in the two phase region will depend dramatically on the system size (calculations with different number of particles become convenient to check the quality of the phase equilibria results) | * The simulation results in the two phase region will depend dramatically on the system size (calculations with different number of particles become convenient to check the quality of the phase equilibria results) | ||
== Direct simulation of the two phase system | |||
== Direct simulation of the two phase system in the [[canonical ensemble]] == | |||
== Gibbs ensemble Monte Carlo for one component systems== | == Gibbs ensemble Monte Carlo for one component systems== | ||
The [[Gibbs ensemble Monte Carlo]] method is often considered as a smart variation of the standard canonical ensemble procedure (See | The [[Gibbs ensemble Monte Carlo]] method is often considered as a 'smart' variation of the standard canonical ensemble procedure (See Ref. 1). | ||
The simulation is, therefore, carried out at constant volume, temperature and number of particles. | The simulation is, therefore, carried out at constant volume, temperature and number of particles. | ||
The whole system is divided into two non-interacting parts, each one has its own simulation | The whole system is divided into two non-interacting parts, each one has its own simulation | ||
box with its own [[periodic boundary conditions]]. | box with its own [[periodic boundary conditions]]. | ||
This separation of the two phases into different boxes is in order to suppress any influence due to | This separation of the two phases into different boxes is in order to suppress any influence due to interfacial effects. | ||
The two subsystems can interchange volume and particles. The rules for these interchanges are | The two subsystems can interchange volume and particles. The rules for these interchanges are | ||
built up so as to guarantee conditions of both chemical and mechanical equilibrium between | built up so as to guarantee conditions of both chemical and mechanical equilibrium between | ||
Line 133: | Line 163: | ||
== Mixtures == | == Mixtures == | ||
=== Symmetric mixtures === | === Symmetric mixtures === | ||
== See also== | == See also== | ||
*[[Gibbs-Duhem integration]] | *[[Gibbs-Duhem integration]] | ||
==References== | ==References== | ||
#[http://dx.doi.org/10.1080/00268978700101491 Athanassios Panagiotopoulos "Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble", Molecular Physics '''61''' pp. 813-826 (1987)] | |||
[[category: computer simulation techniques]] | [[category: computer simulation techniques]] |