Computation of phase equilibria

Thermodynamic equilibrium implies, for two phases $\alpha$ and $\beta$:

• Equal temperature: $T_{\alpha }=T_{\beta }$

• Equal pressure: $p_{\alpha }=p_{\beta }$

• Equal chemical potential: $\mu _{\alpha }=\mu _{\beta }$

The computation of phase equilibria using computer simulation can follow a number of different strategies.

Independent simulations for each phase at fixed $T$ in the canonical ensemble

Simulations can be carried out either using Monte Carlo or Molecular dynamics techniques. 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

- 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 moderate to high density region (where the liquid phase should appear)

- In these simulations we can compute for each density (at fixed temperature) the values of the pressure and the chemical potentials (for instance using the Widom test-particle method)

A quick (and dirty?) method

Using the previously obtained results the following somewhat unsophisticated procedure can be used to obtain a first inspection of the possible phase equilibrium.

Fit the simulation results for each branch by using appropriate functional forms:

$\left.\mu _{v}(\rho )\right.;p_{v}(\rho );\mu _{l}(\rho );p_{l}(\rho )$

Use the fits to build, for each phase, a table with three entries: $\rho ,p,\mu$, then plot for both tables $\mu$ as a function of $p$ and check if the two lines intersect. The crossing point provides (to within statistical uncertainty, the errors due to finite size effects, etc.) the coexistence conditions.

Improving the dirty method

It can be useful to take into account classical thermodynamics to improve the previous analysis. This can be useful because is is not unusual have large uncertainties in the results of the properties. The basic idea is to use thermodynamic consistency requirements to improve the analysis.

Methodology in the NpT ensemble

For temperatures well below the critical point, it could be useful (when the calculation of the chemical potential of the liquid phase Widom test-particle method gives good precision) the following strategy.

• Perform $NpT$ simulation of the liquid phase taking the pressuresas $p\simeq 0$

• Compute the chemical potential of the liquid phase $\mu _{l}(N,T,p=0)$

Van der Waals loops, in the canonical ensemble

Gibbs ensemble Monte Carlo

This method is often considered as a smart variation of the standard canonical ensemble procedure. The simulation is, therefore, carried out at constant volume, temperature and number of particles. Hosever the whole system is divided in two non-interacting parts (each one has its own simulation box with periodic boundary conditions).

The basic idea is to separate the two phases in different boxes to supress the interfacial effects.

The two subsystems can interchange volume and particles. The rules of these interchanges are built up to guarantee that the conditions of chemical and mechanical equilibrium betwwen the two phases are satisfied.

References

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