Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line

Abstract

Thermodynamic integration along a path that coincides with the saturation line is proposed as an efficient means for evaluation of phase equilibria by molecular simulation. The technique allows coexistence to be determined by just one simulation, without ever attempting or performing particle insertions. Prior knowledge of one coexistence point is required to start the procedure. Integration then advances from this state according to the Clapeyron formula—a first-order ordinary differential equation that prescribes how the pressure must change with temperature to maintain coexistence. The method is unusual in the context of thermodynamic integration in that the path is not known at the outset of the process; results from each simulation determine the course that the integration subsequently takes. Predictor–corrector methods among standard numerical techniques are shown to be particularly well suited for this type of integration. A typical integration step along the saturation line proceeds as follows: An increment in the temperature is chosen, and the saturation pressure at the new temperature is ‘‘predicted’’ from previous data (the initial coexistence datum and/or previous simulations). Simultaneous but independent NPT simulations of the coexisting phases are initiated at the said conditions. Averages taken throughout the simulations are repeatedly used to ‘‘correct’’ the estimate of the pressure to convergence. Thus strictly the pressure is not fixed during the simulation. Vapor–liquid coexistence of the van der Waals model is first used to study the numerical integration method without the complications of molecular simulation. In a second application the phase envelope of the Lennard-Jones model fluid is computed, and many variations of the technique are examined. Overall, the results are remarkably consistent and in agreement with previous simulation studies. Difficulty is encountered upon approach of the critical point, but, by artificially coupling the simulation volumes, the method remains effective in this regime so long as a suitably small integration step is employed. Many extensions and improvements of the technique are discussed.