J. Phys. Chem. B 2006, 110, 19601-19604
19601
Competing Evaporation and Condensation Processes during the Boiling of Methane Dirk Zahn* Max-Planck Institut fu¨r Chemische Physik fester Stoffe, No¨thnitzer Strasse 40, 01187 Dresden, Germany ReceiVed: April 10, 2006; In Final Form: July 31, 2006
The atomistic mechanism of the boiling of methane is explored from molecular dynamics simulations. The liquid f vapor transition is initiated by local density fluctuations resulting in a nanometer-sized domain that exhibits both liquid and vapor characteristics. Though the rates of evaporation and condensation events increase dramatically in this area, the overall balance exhibits only a marginal net rate of evaporation. Growth of the precritical domain leads to the nucleation of a vapor phase in which isolated methane molecules are confined by a liquid-vapor interface. After crossing the transition state, the system experiences progressive destabilization of the liquid phase and the evaporation processes clearly outnumber the condensation events.
Introduction Liquid-vapor phase transitions are phenomena of fundamental importance both in nature and for technical applications. Despite their great relevance, our understanding of the underlying mechanisms is still quite limited. This particularly applies to the very early stage of the transformation (i.e., phase nucleation) and subsequent growth. The investigation of such processes requires insights at the atomistic level of detail, which is hard to access from experiment. A new perspective for exploring liquid-vapor transitions was recently demonstrated from molecular dynamics simulations of the boiling of water.1 Therein employing the transition path sampling technique2,3 allowed investigations going far beyond classical nucleation theory. While our recent study was related to the boiling of water, in the present work we shall focus on the liquid-vapor transition of methane. Though methane molecules are of approximately the same mass and size as water molecules, methane lacks the ability to form hydrogen bonds and hence does not display the anomalies observed for water. Indeed, below 111.6 K methane can be considered as a prototype for a “normal” liquid. Moreover, liquid-vapor transitions of methane and alkanes in general are of immense technical importance. While the abundant literature of experimental studies is focused on macroscopic investigations and distillation process optimization, molecular aspects were addressed in a series of theoretical work. The latter field was pioneered by Siepmann, who developed specific model parameters for describing the liquid and vapor phases of alkanes within atomistic simulations.4 This allowed the computation of liquid-vapor phase diagrams and critical temperatures in good agreement with experimental findings.5,6 In addition, Chen et al. explored the nucleation of liquid alkanes from the vapor, applying the commonly used approach of investigating aggregate free energy levels as a function of size.7 While such free energy methods offer important insights concerning static properties, a direct investigation of dynamic features such as competing condensation and evaporation events remains elusive. This limitation may be overcome1 by employing the transition path sampling (TPS) scheme.2,3 Thus far, the TPS approach was used for investigating a variety of phase Corresponding author. Telephone: +49 (0) 351 4646 4205. Fax: +49 (0) 351 4646 3002. E-mail:
[email protected].
transitions.1,8-13 All of these studies benefited from not requiring predefinition of a putative reaction coordinate. As a consequence, the mechanistic analysis may be based on a manifold of unbiased dynamical pathways, each reflecting a phase transition event. The only prerequisite to TPS is an initial trajectory of the transformation used as the starting point for iterative sampling of further pathways. The first trajectory does not need to be a likely one, as the TPS iterations represent a Monte Carlo type sampling of trajectory space and will hence evolve toward the favored regime of transition routes. Simulation Details Our simulation system is composed of 2000 methane molecules described by the all-atom force field of Siepmann and Martin.4 Therein each methane molecule is mimicked by a single effective particle that experiences van der Waals interactions only. Periodic boundary conditions are applied to mimic an infinite liquid (vapor) phase. In the course of the vaporization process, the size of the simulation cell increases dramatically, which motivated the choice of the box volume as an order parameter for evaluating the progress of the transition during the TPS iterations. For the molecular dynamics simulations, a relatively small time step of 0.5 fs was found to be appropriate to ensure time-reversibility, which is of crucial importance for performing TPS simulations. In analogy to our recent study of the boiling of water,1 the initial liquid-vapor transition pathway of methane was prepared from a high-temperature (500 K above the boiling point) vaporization trajectory. From this, snapshots were taken and propagated in both directions of time applying 111.5 K and 1 atm as constant temperature and pressure, respectively. These conditions reflect the liquid-vapor phase coexistence of the used all-atom model of methane4 and deviate from the experiment by only 0.1 K. It is of vital importance to ensure that the liquid-vapor transition is not affected by volume fluctuations related to the constant pressure (and temperature) algorithm used in the molecular dynamics simulations.14 Before the production runs, various barostat relaxation times had been tried: τ ) 1, 5, 10, 25, and 50 ps. Using 1 ps, the barostat causes extreme fluctuations of the simulation cell and the system tends to oscillate from liquid to vapor and vice versa. This reflects a too-strong coupling of the simulation cell volume to the constant
10.1021/jp062221i CCC: $33.50 © 2006 American Chemical Society Published on Web 09/02/2006
19602 J. Phys. Chem. B, Vol. 110, No. 39, 2006
Zahn pressure algorithm. For the larger relaxation times investigated, this phenomenon is not observed. However, using τ ) 50 ps, we found the liquid f vapor transition to be slowed, indicating a choice of a too-high barostat relaxation time. Within the range 5 ps e τ e 25 ps, the transition pathways were qualitatively equivalent. This eventually motivated the choice of τ ) 10 ps for the production runs. Results
Figure 1. Snapshots taken from a representative trajectory illustrating the boiling of methane. Red corresponds to the bulk liquid; methane molecules having no neighbors closer than 0.75 and 1.0 nm are highlighted in green and blue, respectively. (a) Liquid state, (b, c) early stage of vapor nucleation, and (d) formation of the transition state. The latter is compared to (e) the configuration of maximal difference between vaporization and condensation events. (f) Snapshot from methane in the vapor state. Please note that the two-dimensional projection of the pictures makes the vapor domains (b)(e) appear somewhat smaller than they actually are. Additional representation of the transition state configuration (e) is illustrated by the dotted yellow curve.
A typical signature of deriving the initial transition pathway from high-temperature trajectories is related to the strong thermodynamic driving of the process. The latter causes the liquid f vapor transformation to occur in a highly collective, explosive manner. The TPS iterations are hence started in a regime of trajectory space related to spinodal decomposition. However, within about 25 iterations, we observed pathway evolution to the preferred transition mechanism involving phase nucleation and growth as described in the following paragraphs. A characteristic boiling trajectory taken out of 400 transition pathways harvested after convergence to the favored mechanism is illustrated in Figure 1. Therein a color code is applied to mark methane molecules of different coordination. In the liquid state, the first coordination sphere of each methane molecule consists of about 10 molecules at an average C-C distance of roughly 0.5 nm. While the red balls indicate the center-of-mass of “liquid” methane molecules, green and blue coloring is used to highlight methane molecules having no neighbors closer than 0.75 and 1.0 nm, respectively. A snapshot of liquid methane is shown in Figure 1a (t ) -50 ps). In liquid methane at 111.5 K, separation of single molecules by 0.75 nm from the solvent occurs only rarely. Separation by 1.0 nm requires the formation of a vapor domain and is not observed at all in the liquid state. The boiling process is initiated by local density fluctuations resulting in a domain of somewhat lower density compared to the bulk liquid (Figure 1b, time set to zero). In this region, no persistent vacuum regions could be observed. Emerging cavities are instead almost immediately filled by adjacent methane molecules, whose translocation leaves new voids behind. This picture is quite different from the vapor nucleation observed for the boiling of water. In water, the network of hydrogen bonds may allow clathrate-like structure motifs in which cavities are temporarily (around 1 ps) stabilized.1 Methane does not form hydrogen bonds and hence lacks the ability to form such cage structures. Monitoring the formation (f) and decay (r) of coordination shell separation by more than 0.75 nm, we observed a large increase of both kf0.75 forward and kr0.75 backward rates (Figure 2). While the balance kf0.75 - kr0.75 is positive, most of the net production of under-coordinated methane molecules is canceled by competing reassociation to the liquid. However, at one point the system crosses a critical level of destabilization of the liquid state, causing a clear dominance of evaporation. This outnumbering of condensation events is also reflected by the different rates kf1.0 and kr1.0 corresponding to neighbor separation by 1.0 nm or more. Indeed, the related curves are qualitatively very similar, which we take as an indicator that the somewhat arbitrary choice of the distance delimiters does not affect the analysis of vaporization versus condensation. A selection of boiling pathways was scanned for transition states (Figure 3). Following the definition of Du et al.,15 the transition state ensemble is given by all configurations in real space which, after assigning randomly generated, Boltzmann distributed velocities, evolve toward the liquid and the vapor state at equal probability. The boiling trajectories were found
Competing Processes during the Boiling of Methane
Figure 2. Rates k0.75 and k1.0 as a function of time computed for the same boiling trajectory as illustrated in Figure 1. Forward and backward rates corresponding to evaporation (f) and condensation (r) are shown as solid and dotted curves, respectively. The gray lines indicate the times at which the snapshots (b), (d) and (e) in Figure 1 were taken. All rates were normalized with respect to the total number of molecules in the simulation cell. In the stable liquid and vapor states, the forward and backward rates are (on average) equal. The maximum deviation of the two competing rates, denoted as snapshot (e), does not correspond to the transition state (TS), which is instead observed at an earlier stage of the boiling process.
to cross the transition state regime only once. For the vaporization pathway illustrated in Figure 1, the transition state configuration was identified at t ) 7.25 ps (Figure 1d). From the rates plotted in Figure 2, the transition state may be seen to be located at the inset of the divergence of evaporation and condensation. At first sight of Figure 2, one might have expected the transition state to be located at the point of maximal difference of evaporation and condensation rates (Figure 1e). However, the apparent symmetry of kf - kr does not account for the transition state. The snapshots of maximal kf - kr reflect loosely bound agglomerates of still compact but no longer liquid methane molecules. Indeed, the transition state configuration
J. Phys. Chem. B, Vol. 110, No. 39, 2006 19603 Figure 1d corresponds to a critical level of destabilization of the liquid state, while snapshot Figure 1f reflects the result of post-critical destabilization (i.e., the decomposition and evaporation of the condensed phase). Further transition state configurations obtained from other boiling trajectories are shown in Figure 3. Exploring the manifold of possible transition states helps demonstrate the convergence of the TPS procedure. For this purpose, we compared the transition state related to the initial trajectory to the transition states identified in the pathways obtained after a series of TPS iterations. After 25 iterations, the collective destabilization of the liquid observed in the first trajectory evolved toward the expected heterogeneous picture of a localized vapor nucleus embedded in the liquid. Apart from this convergence to the favored transition mechanism, the variation of the shapes and positions of the nuclei in the course of further TPS iterations demonstrates a reasonable sampling of the relevant transition state regime and like this a proper account of the overall manifold of boiling trajectories. In the course of further time propagation, the system is fully transferred to the vapor state (Figure 1f) within 40-60 ps. As the phase transition is investigated at the liquid-vapor coexistence conditions, it is in principle possible to reverse our trajectories and explore the condensation process by reading Figure 1 backward from panel f to panel a. Unfortunately, this task is complicated by technical aspects related to the constant pressure algorithm. The latter adjusts the volume of the simulation cell to maintain the desired pressure. As the number of particles is rather small compared to macroscopic systems, instantaneous averages are not very accurate and both the pressure and the cell volume are subject to fluctuations. In the liquid state, the system’s compressibility is quite low and so are the volume fluctuations. However, in the vapor state the compressibility is dramatically larger, causing much heavier volume changes. While this issue does not affect the time reversibility of our simulations, we feel that the box fluctuation should nevertheless be considered as an artifact with respect to
Figure 3. Series of transition states computed for methane boiling trajectories obtained from the TPS simulations. The trajectory of iteration #1 crosses an unfavorable intermediate derived from a spinodal decomposition pathway at high temperature. This feature disappears during trajectory evolution to the favored mechanism of boiling at 111.5 K. The pathways produced in later TPS iterations clearly show the preference of vaporization from nucleation and growth. The trajectory illustrated in Figures 1 and 2 corresponds to iteration #400.
19604 J. Phys. Chem. B, Vol. 110, No. 39, 2006
Zahn
Figure 4. Transition state depicted from Figure 1e and illustrated for a periodic boundary scenario as used in our simulations. While classical nucleation theory predicts a single critical nucleus of infinite radius, our model results in an infinite number of finite nuclei.
a macroscopic system. As a consequence, the underlying mechanisms of the initial stage of vapor f liquid nucleation as observed from our simulations are of questionable transferability to the real world. While the sketches of the liquid T vapor trajectories related to the first steps of methane aggregation from the vapor might be affected by such artifacts, it is nevertheless worthwhile to discuss later stages of the condensation process from reversing the boiling trajectories. Along this line, condensation may be seen to involve the formation of agglomerates (Figure 1f) whose fate is determined by the competition of condensation and vaporization events. To survive, the loosely bound aggregates shown in Figure 1f need to organize into a more compact liquid domain as indicated in Figure 1e. From a vapor f liquid perspective, crossing the transition state hence implies both growth of the agglomerate of condensed methane and selforganization into a liquidlike form of larger density. Comparing Figures 1e,f and 2, the rearrangement of the agglomerates into more compact structures may be suggested to account for the drop in the vaporization rate and the net increase in condensation (Please note that the curves of forward kf and backward kr rates in Figure 2 must be exchanged upon time inversion). It is worthwhile to discuss our transition path sampling molecular dynamics simulations of the liquid-vapor phase transition of methane with respect to classical nucleation theory. In classical nucleation theory, the phase transition is assumed to occur via a single nucleus, typically taken to be spherical. In this picture, the transition state may be associated to a critical nucleus of radius rcrit. The latter may be calculated from simplified energy considerations of (a) an unfavorable energy term proportional to the surface of the nucleus originating from the formation of a phase interface and (b) a volume term given by the differences in chemical potential of both phases. At liquid-vapor coexistence, the second energy term vanishes and classical nucleation theory leads to a critical nucleus of infinite size and infinite energy of formation. This is in clear contradiction to the observation of the liquidvapor transition in our molecular dynamics simulations. To explain this discrepancy, it is educative to illustrate the effect of using finite simulation models and periodic boundary conditions. For this purpose, we replot the transition state snapshot of Figure 1e including a series of periodic images in Figure 4. Rather than a single critical nucleus of infinite size, our simulations lead to an infinite number of finite phase nuclei. Taken separately, each of the latter of course corresponds to a finite energy of formation. The observation of such transition state configurations of both finite size and energy barrier must
be considered as an effect of the finite size of the simulation box. Again, this limits the transferability of our simulations related to the later stages of the methane boiling process to the real world. Conclusions In conclusion, we presented a molecular dynamics study of the boiling of methane that revealed detailed insights in the nucleation and growth of vapor domains. This process is initiated from local density fluctuations and is accompanied by a large number of both vaporization and condensation events. The competition of the two processes remains almost balanced until the system crosses the transition state regime. At post-critical destabilization of the condensed phase, vaporization processes clearly outnumber condensation events and the entire simulation cell is transferred to the vapor state. The boiling mechanism of methane differs significantly from the liquid-vapor transition of water. The latter was observed to involve small cavities, which are stabilized by hydrogen bonding “around” the vapor bubble.1 Such clathrate-like structures do not exist in liquids which are dominated by van der Waals interactions. While the boiling mechanisms of polar compounds may be expected to depend on the specific molecular structure and charge distribution, the boiling mechanism of methane should be transferable to other liquids which predominantly interact via van der Waals forces. References and Notes (1) Zahn, D. Phys. ReV. Lett. 2004, 93, 227801. (2) Bolhuis, P. G.; Dellago, C.; Chandler, D. Faraday Discuss. 1998, 110, 421. (3) Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. J. Chem. Phys. 1998, 108, 1964. (4) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569. (5) Siepmann, J. I.; Karaborni, S.; Smit, B. Nature 1993, 365, 330. (6) Smit, B.; Karaborni, S.; Siepmann, I. J. J. Chem. Phys. 1994, 102, 2126. (7) Chen, B.; Siepmann, J. I.; Oh, K. J.; Klein, M. L. J. Chem. Phys. 2002, 116, 4317. (8) Pan, A.; Chandler, D. J. Phys. Chem. B 2004, 108, 19681. (9) Zahn, D.; Leoni, S. Phys. ReV. Lett. 2004, 92, 250201. (10) Zahn, D. J. Solid State Chem. 2004, 177, 3590. (11) Moroni, D.; ten Wolde, P. R.; Bolhuis, P. G. Phys. ReV. Lett. 2005, 94, 235703. (12) Zahn, D.; Grin, Y.; Leoni, S. Phys. ReV. B 2005, 72, 064110. (13) Zahn, D.; Hochrein, O.; Leoni, S. Phys. ReV. B 2005, 72, 094106. (14) Melchionna, S.; Cicotti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533. (15) Du, R.; Pande, V. S.; Grosberg, A. Y.; Tanaka, T.; Shakhnovich, E. S. J. Chem. Phys. 1998, 108, 334.
