The Relation of Interface Properties and Bulk Phase Stability

4688

J. Phys. Chem. B 2009, 113, 4688–4697

The Relation of Interface Properties and Bulk Phase Stability: Molecular Dynamics Simulations of Carbon Dioxide T. Kraska,*,† F. Ro¨mer,† and A. R. Imre‡ Institute for Physical Chemistry, UniVersity Cologne, Luxemburger Str.116, D-50939 Ko¨ln, Germany, and Simulator Laboratory, KFKI Atomic Energy Research Institute H-1525 Budapest, POB 49, Hungary ReceiVed: October 5, 2008; ReVised Manuscript ReceiVed: January 20, 2009

The limit of metastability, the so-called spinodal, is calculated for pure carbon dioxide by molecular dynamics simulation. The determination of the spinodal is based on properties of the liquid vapor interface using a recently developed method. This method relates the tangential pressure component through the vapor-liquid interface to the van der Waals loop in the two-phase region of the phase diagram. By application of the thermodynamic stability criteria, the location of the spinodal can be determined. The spinodal determined in this way is called interface spinodal here. Furthermore, the simulation provides equation of state properties in the complete metastable region of the phase diagram. The performance of different correlation equations for the density and the pressure tensor profiles with respect to the estimation of the spinodal is compared. It has been found that the interface spinodal coincides with the thermodynamic mean field spinodal within some reasonable deviation. Finally the influence of the size of the simulation box on the spinodal properties is investigated showing that the temperature-density spinodal data are independent of the interface thickness. Additional simulations using a Lennard-Jones fluid confirm these results over a range of 1.5 orders of magnitude for the systems size. A further result is that interface systems require a very long simulation time in order to obtain reliable results. Introduction Recently it has been shown that the continuously changing properties through an interface can be used to predict the fluid-fluid stability boundaries of homogeneous three-dimensional phases.1 The link between the interface and the homogeneous system is accomplished by a transformation of the data set given by the tangential pressure tensor component and the corresponding density at the same position in the interface. The transformed function is similar to the so-called van der Waals loop through the metastable and unstable region of the phase diagram. This link between interface properties and isotherms of homogeneous phases allows the determination of several properties in the two-phase region from the simulation of an interface. One can, for example, directly determine the thermodynamic stability limit also called spinodal from the interface. Here we call this spinodal “interface spinodal” because the term “surface spiondal” is already used. Technically, the spinodal can be handled as the very final limit of overheating and expansion of a liquid. Before reaching the spinodal other stability limits can appear such as kinetic stability limits, heterogeneous nucleation limits in the case of impurities, or spinodal nucleation, where the activation barrier is in the order of the thermal fluctuations.2 The regime in which classical nucleation theory3-5 is valid is actually a narrow strip close to the binodal in the coexistence region. Close to the spinodal curve a small strip of so-called spinodal nucleation exists. While within classical nucleation the critical nucleus is compact with a sharp interface, in the regime of spinodal nucleation it is rather a ramified cluster without sharp liquid-vapor but diffusive interface2 that also * To whom correspondence should be addressed. E-mail: t.kraska@ uni-koeln.de. † University Cologne. ‡ KFKI Atomic Energy Research Institute.

has consequences for the definition of a cluster in molecular dynamics simulation. In the classical picture based on a mean field equation of state the fluctuations diverge at the spinodal leading to an infinite critical cluster size because the interface of the cluster also diverges toward the spinodal. This divergence becomes dominant only close to the spinodal while in some distance the critical cluster size decreases with increasing supersaturation. In classical nucleation theories the spinodal is not present, and hence the critical nucleus always shrinks when the supersaturation rises. In any case close to the spinodal, the already highly metastable liquid will loose the stability completely, and therefore the liquid experiences sudden and fast pressure change leading to an explosion-like boiling also called flash boiling.6 Spinodal data can therefore be employed for the worst case scenario because in the vicinity of the spinodal explosive boiling is most distinctive. Here we present new results for carbon dioxide, which is a widely applied supercritical fluid used in extraction and particle production processes.7 Such processes typically take place far away from the equilibrium state which requires the knowledge of the stability limits and the description of the states in the metastable region. Recent plans to inject the combustion gas carbon dioxide subsurface, where it can be retained for geological periods of time, would also produce huge amount of high pressure and high temperature carbon dioxide, where safety calculations concerning possible gas leakage require the knowledge of the stability limits.8 Methods Potential Model. For the simulations we employ the rigid EPM2 potential model for carbon dioxide.9 Within this model the carbon dioxide molecule is represented by three LennardJones (LJ) sites and three partial charges at the same positions.

10.1021/jp808789p CCC: $40.75  2009 American Chemical Society Published on Web 03/10/2009

MD Simulations of Carbon Dioxide

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4689

TABLE 1: Potential Parameters of the EPM2 Carbon Dioxide Model qC/e qO/e σO/nm σC/nm O/K C/K dCO/nm dm/nm shake tolerance

0.6512 -qC/2 0.3033 0.2757 80.507 28.129 0.1149 0.1980 10-3

To keep the molecule linear and rigid we perform force shifting according to Ciccotti et al.10 and map the masses of the molecule into two mass points which are kept nearly in constant distance by the Shake algorithm.11 The choice of this model is two-fold: first we need a model with sufficient accuracy to describe supercritical fluids and second we need a model compatible to that of large organic molecules in the continuing work. For the interaction between two LJ sites the Lorentz-Berthelot mixing rules are used. For the long-range Coulomb interactions the Ewald sum is employed. For all ensembles except of NpT the particle mesh Ewald (PME) method with a grid distance of 0.1 nm and with fourth-order interpolation is utilized to enhance the simulation speed. Using PME with these parameters leads to very good results for the surface tension for a given potential model.12 The large cutoff of the real-space contribution of the Ewald sum, which corresponds to the LJ cutoff, leads to a relatively small long-range correction.13 The LJ interactions are truncated at a cutoff distance of 1.7 nm, which is 5.6-6.2 times the LJ diameter of the sites and usually larger than the interface thicknesses obtained in our calculations (see below). This large cutoff radius has been chosen in order to minimize cutoff-effects of properties for which no long-range corrections are available. Simulation. All simulation presented here are performed with the MD code MOSCITO.14 We use a time step of 1 fs and write out every 0.2 ps the system configuration. Setting up a film simulation requires several steps in order to keep the film stable during equilibration. For better efficiency of the equilibration simulation it is useful to employ a sequence of different ensembles. To simulate liquid films in equilibrium with a vapor phase we start with a simulation box at a density slightly below the corresponding liquid density at the given temperature and perform an NVT run lasting 1 ns for equilibration. To keep the temperature constant in the equilibration run the Berendsen thermostat15 is used. This run is followed by a NpT run for 2 ns using the Berendsen barostat for given equilibrium vapor pressure and temperature in order to stabilize the system. Finally, the simulation box is elongated three times the box length in one direction for which we chose the z coordinate. For the investigation of the influence of the box length on the properties elongations (Lz) up to nine times the original box length are used. The system is then again equilibrated in an NVT run for 2 ns. After equilibrium is reached, i.e., vapor and liquid density are constant and temperature of both phases are equal, we finally run a constant energy or NVE simulation for 1 ns. We investigate two system sizes, a small system containing 1000 and a big system containing 2197 CO2 molecules. For the small system the initial box size is 3.91 nm × 3.91 nm × 11.7 nm and 6.04 nm × 6.04 nm × 18.12 nm for the big system. These dimensions correspond to the above-mentioned elongation of Lz ) 3. Data Analysis. Once we have obtained the configurations of the stable liquid film in equilibrium with its vapor phase for a sufficient period of time they are analyzed with respect to the density and pressure tensor profiles. The density profiles are

fitted with two different functions, namely, the tangent hyperbolic16 and the error function profile17 given by

F(z) ) 0.5(Fl + Fv) - 0.5(Fl - Fv) tanh

( 2(z d- l) ) (1a)

and

F(z) ) 0.5(Fl + Fv) - 0.5(Fl - Fv) erf

(

√π(z - l) d

)

(1b)

respectively. Here Fl and Fv are the equilibrium liquid and vapor densities, l is the midposition of the interface, and d is the thickness parameter. The forms of the two equations are similar, but it should be noted that besides the different functions (tanh and erf) the constant in the argument is different. The tanh profile is based on the density gradient approximation and gives good results for mean field simulations such as the Lattice-Boltzmann simulation,18 while the erf profile is a result of the capillary wave approach.17 This capillary wave approach exists19-21 besides the assumption of an intrinsic interface between the liquid and the vapor phase. In this approach it is assumed that the interface is sharp as the Gibbs dividing interface and thermal fluctuations lead to a wider interface density profile. This idea has also been applied to already diffusive interfaces where the capillary waves lead to a further broadening.22 As mentioned above, one result of the capillary wave approach is the error function density profile (eq 1b), which has been obtained by averaging over the surface partition function.17 In the context of this work here capillary waves may have an effect on the spinodal data. It has been found23 that in several cases such as very thin films, the presence of real walls, or external fields or the use of the grand canonical ensemble the interface width can significantly depend on the system size. However these cases are not relevant here, still the remaining influence of the system size may exist and is discussed below. We compare both profiles (eqs 1a and 1b) in order to analyze its influence on the interface spinodal estimation. Numerical analysis shows that the d value of the tanh density profile function represents a 12-88 interface. This means that the thickness is defined as the distance between the z coordinates at which the density is Fv + 0.12(Fl - Fv) and Fv +0.88(Fl Fv). It should be noted here that experimental data are usually obtained for the 10-90 surface. Further analysis of eq 1a shows that the 10-90 interface is approximately 10% thicker than the 12-88 surface as obtained by the parameter d of the tanh function. This is a non-negligible effect and should be included in the analysis of the data especially when comparing different data. The erf density profile function (eq 1b) gives the interface thickness between 10.5 and 89.5% of the density difference, which is closer to the 10-90 interface. The pressure tensor is calculated using the Iring-Kirkwood convention.24 Its components can be divided into two contributions: the kinetic part is the same for the normal and the tangential component, and the configuration part which is different for the normal and the tangential component. The configuration part is calculated by the usual virial expressions as described in a previous paper.1 The kinetic part of the pressure component can be correlated to the same type of function as the density profile. The normal pressure tensor component pN

4690 J. Phys. Chem. B, Vol. 113, No. 14, 2009

Kraska et al.

perpendicular to the interface is constant through the interface which can be shown by elasticity theory.25 It is furthermore equal to the equilibrium vapor pressure. We fit the difference ∆NTp ≡ pN - pT to a function given by1

∆NTp )

(a - b)(1 + (tanh(Z))2) - 2(a + b) tanh(Z) (cosh(Z))2

(2) with Z ) 2(z - l)/d and the fit parameters a and b. The form of this function is the result of the application of the equation of Fuchs26

∆NTp ) -q(F(z)′′F(z) - (F(z))2)

(

2π∆F(z - l) exp (-y2) (Fl + Fv 2d3 m∆F2 (exp(-y2))2 ∆F erf (y)) d2

)

(4)

with ∆F ) Fl - Fv and y ) (π)1/2 (z - l)/d. Here the parameter m is introduced in order to get a better correlation of the simulation data. The exact expression for the application of eq 3 to eq 1b is obtained by setting m ) 1. The property of interest here, the interface spinodal, is calculated from the tangential pressure component. This approach is based on the idea that the tangential pressure component plotted vs the density at the same position in the interface gives a function which is related to the van der Waals loop of an equation of state in the two-phase region. A first hint that a density exists in the interface which corresponds to nonstable states can already be found in the paper on capillarity by van der Waals.27 It was also discussed by Bakker who drew a van der Waals like loop of the pressure-density relation in the interface.25 However the discussion by Bakker was not correct because he oversaw the small maximum on the tangential pressure on the vapor side and he assumes a discontinuity of the isotherm at the binodal. Recently it has been shown that the tangential pressure component can actually be transformed into a van der Waals like loop by plotting it versus the density at the same position (i.e., z coordinate here) in the interface.1 This requires a linear transformation of the tangential pressure

3 p(z) ) pN - (pN - pT(z)) 2

psp,liq ) pvap - s

(3)

to the tanh density profile where q is a parameter that scales the resulting function to the pressure dimension. The value of q can in principle be calculated from an integral over the interaction potential and hence should be substance- and temperature-dependent parameter. The form of eq 2 ensures the correct limit toward the bulk liquid and vapor phase of the tangential pressure component if calculated by pT ≡ pN - ∆NTp where the normal pressure is the equilibrium vapor pressure. If one applies eq 3 to the erf profile eq 1b one obtains

∆NTp ) -q

critical point both pressure components are equal leading to a + c ) 1 and hence p(z) ) pN - c(pN - pT(z)); second assuming a transformation from the tangential pressure in two dimensions to the three-dimensional system the factor c takes the value 3/2. It should be noted that this value influences the interface spinodal pressure only but it does not affect the temperaturedensity diagram. The interface spinodal can then be determined by applying the mechanical stability criterion from the extreme values of the van der Waals loop. Furthermore, it has been shown that one can estimate the spinodal from experimental equilibrium data being the vapor pressure pvap, the surface tension γ, and the interface thickness d by

(5)

This transformation is based on the assumption of a linear composition of the normal and the tangential pressure given by p(z) ) a pN - c pT(z). The coefficients a and c of the two contributions are obtained from two side conditions: first at the

3γ 2d

(6)

which can be derived1 from eq 5 using the mechanical definition of the interfacial tension γ ) ∫(pN - pT(z))dz and setting pN ) pvap. Here s is a so-called shape factor describing the approximate shape of the pT(z) profile on the liquid side of the interface region. If we assume a rectangular peak as approximation for the real peak, we obtain the value s ) 1; for a triangular peak s would be 2. Here s is neglected, i.e., set to unity. It should be noted that eq 6 is only valid significantly below the critical point because it neglects the small but wide maximum in pT on the vapor side of the interface. This peak becomes negligible in sufficient distance to the critical point. Furthermore, it should also be emphasized that eq 6 contains only measurable quantities and does not have any adjustable parameter. It is therefore a useful tool to estimate the liquid spinodal, which cannot be measured directly. Results Simulations are performed for six temperatures ranging from 220 to 280 K covering two-thirds of the liquid region. The triplepoint temperature for carbon dioxide is 216.6 K, and the critical temperature is 304.2 K. A representative snapshot of a simulation of 2197 molecules at 260 K is shown in Figure 1. At too high temperatures the interface becomes so wide that it does not fit into the simulation box. The density and pressure tensor profiles are calculated from 5000 configurations sampled over 1 ns simulation time. Before calculating the profiles the moment of inertia of the liquid film is calculated from the first moment of the density distribution in the liquid film, and then it is shifted to z ) 0. This allows to average over both surfaces of the film. Figure 2 shows the density profiles for 220 K. The simulation data are smooth and can be very well correlated by the tanh and the erf profile on the given scale (Figure 2a). Small deviations are visible in the region of large curvature. The enlargement of these regions shows that the tanh profile seems to be slightly better on the liquid side (Figure 2b) while the erf profile is significantly better on the vapor side (Figure 2c). Especially in the steep region of the profile on the vapor side the deviation of the tanh profile can have a significant effect on the calculation of the density. The interface thickness d obtained from the correlations changes from 0.9 to 1.8 nm and is plotted in Figure 3 vs the temperature. For comparison, the length of the linear CO2 molecule can be approximated as twice the CdO bond, which is 0.23 nm. The d value obtained from the erf profile is consistently a few percent above that of the tanh profile. The same trend is observed for the systems with 2197 molecules. The reason for this systematic difference is the inherent definition of the interface thickness of the tanh- and



Figure 1. Snapshot of a simulation with 2197 CO2 molecules at 260 K. Green, carbon; red, oxygen.

the erf-function as discussed above. Correlation with a power law gives exponents close to -0.5 for N ) 1000 and somewhat lower values for N ) 2197 due to the larger scattering of the data. The pressure tensor components obtained from the simulations are shown in Figure 4 together with the correlations. The kinetic part of both the tangential and the normal pressure component is plotted in Figure 4a for two temperatures. It gives a profile similar to that of the density profile and can be fitted by the same functions giving the same value for the interface thickness d as the density profile correlation. The difference between the normal and the tangential pressure is plotted in Figure 4b for 220 K. The correlation with eq 2 is quite good; however, as reported earlier,1 the small extreme value on the vapor side of the interface may not be very accurate. This is because the residuals are very small compared to that of the large extreme value on the liquid side. Hence it is necessary in some cases to analyze the vapor side extreme value separately. By use of eq 5 we can calculate the interface spinodals from the extreme values of the tangential pressure tensor. The result is plotted in Figure 5a together with the vapor pressure curve. The open circles are the interface spinodal data for the small system (N ) 1000) obtained from the correlation with eq 2. They agree well with calculations using the Peng-Robinson equation of state.28 The vapor spinodal data at low temperature may be inaccurate, which is related to inaccuracies in the determination of the extreme value in ∆NTp. As just mentioned the small vapor phase extreme value may not be well represented by this function. This effect becomes stronger with decreasing temperature because the liquid peak of the tangential pressure component becomes larger, and hence the vapor peak becomes almost irrelevant in the correlation. The liquid spinodal for the big system is systematically above that of the small system. The reason for this observation is discussed below within the analysis of finite size effects. In Figure 5a we have also added the interface spinodal estimated with eq 6 (marked by a plus sign) using the simulation data for the vapor pressure, the surface tension, and the interface thickness for the small system. These data are consistently a little bit above the spinodal data obtained directly with eq 5. However, the difference is relatively small, and it follows that the value s ) 1 can be regarded as reasonable approximation for the liquid spinodal of carbon dioxide. Therefore we keep this value fixed rather than using it as adjustable parameter. This is consistent with the simulations of argon, which also gave the value s ) 1 over a wide temperature range significantly below the critical point.1 In addition we have applied eq 6 using

experimental data for the vapor pressure, surface tension, and liquid-vapor interface thickness. While the first two properties can be accessed easily, interface thickness has been measured only by Schmidt and Moldover29 very close to the critical point (295.16 K). They measured the correlation length, being half of the interface thickness, by ellipsometry and obtained 0.684 nm corresponding to an interface thickness of approximately 1.37 nm. This gives a spinodal pressure of 4.97 MPa for s ) 1 which is added in Figure 5a. It should be noted in this context again that although this sole point lies close to the calculated spinodal, in general the applicability of eq 6 is weak close to the critical point. The temperature-density plot of the phase diagram is shown in Figure 5b. In this plot there is no systematic difference in the interface spinodal data of the small and the big system. The simulation data are compared to two equations of state, namely, the PR28 and the LK equations of state.30 The latter one has been developed to improve the description of the near critical and the high pressure regions. While both equations of state exhibit good agreement on the vapor-side binodal and spinodal, the PR equation is inaccurate for the liquid binodal. Hence the good agreement of the simulation data and the PR equation for the liquid spinodal in this projection is not meaningful. The LK equation accurately describes the saturated liquid densities and also well represents the liquid spinodal simulation data in the temperature-density projection. This shows that for the simultaneous description of the vapor spinodal and the saturated vapor density a simple equations such as the PR equation is sufficient while for the liquid side more accurate equations of state such as the LK equation are required. On the other hand the prediction of the liquid spinodal pressure by the LK equation (Figure 5a) is not very good. A reason could be that this equation is fitted using equilibrium data especially in the critical region and at high pressure leaving some deviations for the extrapolation in the two phase region at low temperature. We have found similar deviations in the liquid spinodal pressure for other noncubic equations. Such equations can be further improved using the interfacial data of the metastable region as determined here. As discussed above the representation of the vapor peak of the tangential pressure in the interface by the tanh-profile (eq 2) resulting from the application of the Fuchs equation (eq 3) on the tanh density profile (eq 1a) is not always very good, especially at low temperatures. Therefore we here also investigate the same approach, the Fuchs transformation, by using the erf density profile eq 1b which leads to eq 4. In the correlations shown in parts a and b of Figure 6, one can see


Kraska et al.

Figure 3. Interface thickness d obtained from the correlation using the tanh (circles, N ) 1000; triangles, N ) 2197) and the erf profile (squares, N ) 1000) function. The d values of the erf function for N ) 2197 (not shown here) are systematically above the corresponding values obtained from the fit to the tanh-function as for the small system. The elongation is Lz ) 3 for all data points. The curves are correlations with a scaling equation.

Figure 2. Density profile of the carbon dioxide vapor-liquid interface as obtained from MD simulations of the small systems (N ) 1000) (gray) and calculations with the tanh (red, dashed) and the erf (blue, solid) profile functions. (a) Complete density profile. (b) Enlargement of the liquid side. (c) Enlargement of the vapor side.

that this equation represents the tangential pressure profile somewhat better. The parameters of eq 4 including the empirically added parameter m are obtained by a fit to the pN-pT simulation data. The introduction of the parameter m improves the vapor side (Figure 6b). It should be noted that both correlations of the Fuchs-transformed erf function, with and without adjustable m, give values for the parameters, such as d, which are different to the ones obtained form the erf density profile correlation with eq 1b. This basically reflects the quantitative inaccuracy of the Fuchs transformation. In theories for surface tension based on this transformation such as the square-gradient theory of van der Waals,27,31 this inaccuracy is corrected by an adjustable parameter. However, the spinodal obtained from the direct correlation of eq 4 to the pN-pT simulation data including the adjustable parameter m does not show significant deviations to that obtained from eq 3.

With the spinodal we compared so far only one fix-point in the two-phase region to equation of state calculations. However, the determination of the isotherms in the metastable region is also of practical interest. To construct the van der Waals loop we plot in Figure 7 p(z) as function of F(z) using the correlation functions 1a and 2 fitted to the simulation data together with the linear transformation given by eq 5. The resulting data are plotted as solid curves in Figure 7. The isotherms of the LK equation of state are systematically at lower pressure than the simulation data. On the other hand the density scale appears to be the same for the equation of state and the simulation. In addition we plot in Figure 7 also a transformation with eq 5 but change the factor 3/2 to other values in order to get the best possible agreement. The values of this factor c given in Figure 7 vary from 1.8 for 230 K to 2.5 for 270 K. Using a variable factor apparently gives very good agreement. Some deviations on the liquid side are related to general inaccuracies of an equation of state. Since equations of state can only provide extrapolations into the two phase region it remains unclear whether the factor 3/2 in eq 5 is inaccurate or whether the equation of state is inaccurate in the two-phase region. It should be noted that this factor affects the pressure transformation only and not the density scale. Interestingly, fitting the adjusted c parameters to an exponential function of the temperature gives a limiting value for c slightly below 3/2, which may show that this c value is the low temperature limit. As already mentioned eq 6 is valid in sufficient distance to the critical point. Still one hat to keep mind the general inaccuracy of equations of state in the two-phase region. It has been discussed whether van der Waals loops obtained in the Monte Carlo simulation of homogeneous finite size systems are related to the those obtained from an equation of state.36 Therefore, as another test for the data obtained here from the interface simulation one can apply the Maxwell construction. This is a well-known method for the calculation of phase equilibrium for pure substances which is often used as a simple teaching illustration. Of course it can be derived rigorously from the phase equilibrium condition, which is Gm,liq ) Gm,vap for a pure substance. Inserting the relation between the Gibbs energy G and the Helmoltz energy A leads to Am,liq + peqVm,liq ) Am,vap + peqVm,vap with the equilibrium pressure peq which is equal in both phases. Furthermore, using A ) - ∫p dV gives m,vap the underlying expression of the Maxwell construction ∫VVm,liq p(Vm)dVm ) peq(Vm,vap - Vm,liq). Hence there is a relation between the isotherm in the two-phase region and the stable



Figure 4. (a) Profile of the kinetic pressure through the vapor-liquid interface obtained by MD simulations. The curves are correlations with an equation having the structure of eq 1a. (b) Correlation of the tangential pressure profile with eq 2 (curve). The symbols are simulation data.

Figure 5. (a) Pressure-temperature projection of the vapor pressure curve, the vapor spinodal above it, and the liquid spinodal below it. Diamonds, vapor pressure data for N ) 1000; circles:, spinodals for N ) 1000; triangle up, vapor pressure curve for N ) 2197; triangle down, spinodals for N ) 2197; plus sign, liquid spinodal calculated with eq 6 using the MD results for the vapor pressure, surface tension, and interface thickness and s ) 1; hexagon, calculated from experimental data using eq 6 with s ) 1. Curves: calculations with the Peng-Robinson (PR) and the Leonhard-Kraska (LK) equations of state. (b) Coexistence curve and spinodal obtained from simulations for carbon dioxide. Legend as for part a.

Figure 6. (a) Pressure tensor profile correlation with the error function (erf). Solid (red) curve, eq 4 with m ) 1; dash-dotted (blue) curve, eq 4 with adjustable parameter m. (b) Magnification of the vapor-side peak of part a.

phase equilibrium but the course of the isotherm in the twophase region cannot be uniquely determined by the Maxwell construction. Inverting eq 1a and inserting it together with F ) 1/V into eq 2 gives an expression for p(V). After subtracting pN one can calculate the two areas on the liquid side and on the vapor side which should be equal. As a result we find for T e 270 K some differences between the area on the liquid side Oliq and on the vapor side Ovap ranging from (Oliq + Ovap)/ ((| Oliq| + |Ovap|)/2) ) (3-20%, but these deviations appear to be statistical uncertainties because they fluctuate as function of the temperature. This means that once the area on the liquid side and once that on the vapor side is larger. These uncertainties are caused by of the simulation data itself and the inaccuracies due to correlation by eqs 1a and 2. Especially on the liquid side, where the isotherm is very steep, small inaccuracies can lead to such deviations. The variation of the factor c being 3/2

in eq 5 does not affect the Maxwell construction. The same evaluation using earlier simulation data of argon1 gives the same result. Finite Size Effects. It is known that the finite size of a simulation system affects the interface properties such as the interface width and the surface tension. The origin of these effects are capillary waves23 and further oscillating finite size effects in the limit of the very small systems.32,33 These oscillations in the surface tension have been found for LJ and ionic liquid systems with a vapor-liquid interface smaller than (10σ)2. Above that size only nonoscillating finite size effects exists according to the capillary wave approach. The systems investigated here consist of 1000 and 2197 carbon dioxide molecules. The small system has a vapor-liquid interface of (12.9σ(O))2, while the large systems have an interface of (19.9σ(O))2. To analyze the effect of the system size on the


Figure 7. Plot of p(z) as function of F(z) (black, solid curves) calculated by eq 5 in comparison to the isotherms of the LK equation of state (red, dashed curves). The (blue) dot-dashed curves are calculated from the simulation data using eq 5 but with the factor c given in the diagram instead of the factor 3/2.

Figure 8. Interface thickness for different system size as function of the elongation of the simulation box in the z direction.

estimation of the interface spinodal we have performed for both systems sizes simulations with various elongations. In Figure 8 the interface thickness is plotted as function of the elongation Lz. The parameter Lz is the factor of the box expansion relative to the width of the cubic simulation box of initial bulk liquid phase simulation. For the small system with 1000 molecules we find that the interface thickness varies significantly and nonmonotonically with the elongation. Furthermore, no clear trend of the interface width with the elongation is observed as discussed in the literature.23 The big system with 2197 molecules exhibits a much weaker variation of the interface thickness. In Figure 3 the effect of the system size by means of number of molecules and interface area for the same elongation is compared. In most cases, especially at low temperature the interface thickness of the big system is higher than for the small system which is in agreement with the capillary wave approach. For higher temperature, i.e., closer to the critical point, the accuracy of the simulation is worse due to the small elongation Lz ) 3. The capillary wave approach loses its validity when approaching the critical point, i.e., at high temperature. The question is now how these effects influence the estimation of the interface spinodal. To analyze the influence of capillary waves here we focus on the effect of the interface thickness because it is known to be affected by capillary waves. For this purpose a model is required that gives consistent density and pressure tensor profiles. We start with a density profile for carbon dioxide as obtained from the molecular simulation. By application of a Fuchs transformation (eq 3) we calculate the corresponding values for ∆NTp. For the parameter q in eq 3 one can chose one since now we are only interested in the relative trends. From the maxima and minima in ∆NTp(z) we obtain the corresponding z values which are then used to

Kraska et al. calculate the spinodal densities from the density profile. For a systematic study we vary the thickness parameter d in eq 1a. As a result we find that the densities of both the liquid and the vapor spinodal do not depend on the interface thickness. The reason is that if capillary waves lead to a wider interface it widens simultaneously the F(z) and the ∆NTp(z) profiles. In addition, the fact that the correlation of the density profile (Figure 2a) and the kinetic pressure (Figure 4a) give the same interface thickness d also shows that the density profile and the tangential pressure profile widen simultaneously. Therefore the maxima and minima in ∆NTp(z) are at the same density and have the same value, no matter how wide the interface is. As a result the interface spinodal densities in the T(F) diagram are not affected by capillary waves (Figure 5b). To verify these results additional MD simulations of argon films using the LJ potential are performed. We use the argon parameters but with cutoff of 2.5σ (0.85225 nm) in order to be able to calculate much larger systems than it would be possible for carbon dioxide. The short cutoff affects the properties of the system significantly34 which has to be taken into account. The simulations here are for 77 K corresponding to 96 K for the full LJ fluid, i.e. with a large cutoff of 5 or 6.5σ. The system sizes are N ) 1000 (10.722 nm, 3.574 nm), N ) 2744 (15.468 nm, 5.156 nm), N ) 10648 (24.315 nm, 8.105 nm), and N ) 32768 (34.371 nm, 11.457 nm). This covers 1.5 orders of magnitude. Interface profile data are written on disk in 0.1-ns blocks and averaged and analyzed together with all preceding data. With proceeding simulation time the statistics becomes better and the results become more accurate. The results for the liquid spinodal data of these systems are plotted in Figure 9 vs the simulation time. The results of these simulations confirm that the spinodal density is independent of the box size. The liquid spinodal density converges for all system sizes toward almost the same value (Figure 9a). The absolute difference in spinodal density is negligible compared to the inaccuracy of the simulation and the correlation. It is virtually invisible in the full temperature-density phase diagram (Figure 9b). The slight trend of decreasing liquid spinodal density with increasing systems size is by orders of magnitude smaller than that found by Kaski and Binder.37 For this and other reasons discussed below we conclude that the interface spinodal density is not size dependent. On the other hand we find that the liquid spinodal pressure is size-dependent (Figure 9c). With increasing system size the interface thickness calculated from the density profile rises as shown in Figure 9d. This also widens the pT profile (d from pT profile in Figure 9d) leading to a flattening of the pT peaks in the interface and hence to a less negative spinodal pressure. Interestingly the plot of the liquid spinodal pressure versus the γ/d in Figure 9e exhibits a linear dependence. As a consequence eq 6 turns out to be independent of the system size and can hence be applied to experimental macroscopic data. However the results also show that the one has to use a relatively long simulation time for averaging in order to get reliable results and to verify the effect of system size on the spinodal pressure. For CO2 we find the same trend as observed for argon, namely, that the spinodal pressure for the large system is above that of the small system (Figure 5a). The enormous computational effort for simulations using the electrostatic CO2 model limits us to these two system sizes but agreement with the results of the argon calculations is obvious. Discussion It is useful to classify the different results from the literature for simulation in the nonstable (metastable and unstable) region



Figure 9. Plot of the liquid spinodal data at 77 K for LJ argon systems (cutoff 2.5σ) of different sizes. (a) Liquid spinodal density as function of the simulation time. (b) The same as part a but on the scale of the complete phase diagram. (c) Liquid spinodal pressure as function of the simulation time. (d) Interface thickness parameter d as function of the system size. The top data are obtained from the fit of the tangential pressure by eq 2, while the bottom data are obtained from the fit of the density profile to eq 1a. (e) Plot of the liquid spinodal pressure vs surface tension divided by interface thickness obtained from the density profile. The slope of the linear regression (dashed line) is -2.06.

that are used to obtain spinodals and other stability limits. There are investigations of finite size system mostly by Monte Carlo simulations35-38 but also by molecular dynamics simulations39,40 inside the nonstable region. These simulations give a loop resembling a van der Waals loop, but it is not a continuous function of a homogeneous density. Actually the extrema of this loop correspond to the onset of two-phase coexistence at equilibrium conditions for given system size. Hence in the density range between the extrema of this loop the system consists of a droplet in equilibrium with the vapor, i.e., it is a heterogeneous two-phase system. The fact that one obtains a loop is related to the additional surface energy of the droplets. These extrema have been called “effective spinodals,” but actually they are indicating the densities between which the system is heterogeneous. In the thermodynamic limit of an infinite system the branch between the extreme becomes the tie-line while the extrema reach the coexistence curve. However this approach lacks the appearance of the thermodynamic spinodal which is due to the ill-defined metastability within statistical theory as discussed below. It would mean that there

are no metastable and no unstable regions in a system that may be large, even macroscopic. Another approach is to investigate the van der Waals loop of a homogeneous system through the nonstable region. In contrast to the before-mentioned investigations one tries to avoid phase separation by constraints in order to get an equivalent equilibrated system for the investigated metastable state.41 One constraint can be the period of time. A simulation is started at a density between the coexistence densities, i.e., in the nonstable region. The data gathered in the period of time before the system decomposes is used to calculate the properties of a homogeneous state in the two-phase region. The higher the supersaturation the faster the system decomposes, and hence the less data are accessible for averaging. A second constraint can be applied on space. To hinder the system to decompose one can suppress density fluctuations. This can be accomplished by using small simulation systems or by introducing confinement cells.42-44 In this way one can obtain data of homogeneous states beyond the binodal in the metastable region which can be compared to experimental data. One can even go beyond the spinodal into

4696 J. Phys. Chem. B, Vol. 113, No. 14, 2009 the unstable region which is however not accessible experimentally. The resulting pressure-density plot corresponds to the van der Waals loop and its extrema to the thermodynamic mean field spinodal. As discussed by Langer,45 on the basis of the statistical mechanics one can only obtain a Helmholtz-free energy functional by coarsening, i.e., averaging over a system or a part of a system of a certain size. The question is about the size of the system for averaging. Langer gives two limits: the system should not have linear dimensions noticeably larger than the correlation length to ensure local thermal equilibrium, and the system size should be large compared to the volume of the molecule. In case of a vapor-liquid transition Langer points to the problem that for a low density vapor phase the coarsegraining cell size has to be relatively large to fulfill the above conditions. On the other hand a relatively large cell-size also means that fluctuations are not properly suppressed and phase separation might occur. Consequently a small condensing droplet would fit into a single cell which does not allow the coarsening for the liquid phase. It seems to be impossible to obtain a Helmholtz-free for such system; however, this is the case for a two-phase system in equilibrium only but not for a homogeneous system in the nonstable region. Kaski et al.37 calculated a spinodal by Monte Carlo simulation using a similar idea as Langer to study the influence of the coarse-graining size, i.e., the size of cell into which the system is divided. They find a significant dependence of the spinodal on the cell size when it becomes larger than the correlation length. For large cell sizes the spinodal approaches the coexistence curve. Again, this is because the system is heterogeneous between these spinodals and not homogeneous. Therefore that spinodal was called “effective spinodal” as discussed above. It is not the thermodynamic spinodal that one can obtain from the extrema of a van der Waals loop. These different approaches result from different interests. Statistical theory is valid in the thermodynamic limit which is either infinite large system size and/or infinite observation time since due to the ergodicity, ensemble and time averages are identical in that limit. From experiments it is known that metastable states can exists for sufficient time to gather reproducible data but compared to the thermodynamic limit, namely infinite time, the period of time in which a metastable state exists is negligible. Hence its contribution to the partition function is also negligible. Therefore in statistical theory metastable states do not exist, i.e., the “true free energy”, which was meant according to statistical theory, exhibits no metastable state.46 For not too high supersaturation the existence of a metatstable state, also called time lag or waiting time for the first nucleus, corresponds to the reciprocal nucleation rate.47 Molecular dynamics simulations confirm that even at high supersaturation nucleation rates obtained from the a waiting time for the first nucleus fit into the order of magnitude of nucleation rates obtained from the steady state statistics.48 With increasing supersaturation the waiting time decreases and the nucleation rate increases until the system spontaneously phase-separates when it is the close to the spinodal. This has been shown recently in simulations of LJ argon.48 That work focused on the possibility of nucleation simulations in a NVE ensemble which is why a very high supersaturation was required in order to remove sufficient energy from the system at the initial temperature jump. In some cases one could observe spontaneous phase separation while at lower supersaturation a time lag was observed. This finite time period or time lag is one possible constraint to probe a metastable system as discussed above. The

Kraska et al. other constraint, the introduction of a finite size, can also be accomplished experimentally.49 Actually both constraints are very much the same as introduced in molecular simulations of metastable states toward the mean field spinodal described above, namely, the observation time prior to decomposition of the system and the introduction of space confinement. The question remains what type of spinodal is the interface spinodal calculated here? Since it is obtained from a function that varies continuously with the density, namely, through the interface, it corresponds to free energy functional of a homogeneous system. For each z position in the interface we get one density without significant fluctuation in the xy plane that could cause a phase separation in the xy plane. Therefore it does not correspond to a free energy functional that represents a heterogeneous system between the so-called “effective spinodals” and hence not necessarily exhibit a size dependence. For the spinodal density we give convincing arguments that this is actually the case, by analytical considerations, by simulations for the CO2 system, and especially by simulations of a LJ system investigated over 1.5 orders of magnitude. There is a difference between finite-size effects in bulk systems and in interface systems. In bulk-phase simulations one uses a small piece of a macroscopic system while in interface systems, sufficiently below the critical state, the characteristic length in one dimension, namely the interface thickness, is smaller than the simulation box length. There are still finite size effects in the other dimensions caused by capillary waves; however, these effects cancel for the temperature-density plot as discussed above. So by using two-dimensional interface properties we can at least reduce the size dependent spinodal properties by the spinodal density. The remaining finite size effects on the spinodal pressure can be explained by the onset of capillary waves in the simulation system of given size. Werner et al.50 concluded from simulations of polymer interfaces that the interface width is based on an intrinsic profile related to a mean field theory that is broadened by capillary waves. For small systems the interface spinodal pressure is that of the mean field system in the same way as the intrinsic density profile is that of the mean field system. So for small systems the interface spinodal pressure hence corresponds to the spinodal obtained from an equation of state. With increasing system size capillary waves set in leading to a broadening of the interface and affecting the interface spinodal pressure. However, as indicated in Figure 9e, eq 6 is size independent and can therefore also be used for macroscopic experimental systems. It follows that the interface spinodal does not approach the coexistence curve in the thermodynamic limit as the “effective spinodal” does. Heermann et al.51 have investigated the vapor-liquid nucleation of carbon dioxide modeled by a LJ force field by calculating the compressibility. In mean field approaches the compressibility diverges at the spinodal while it rises to a finite value in non-mean field models. The spinodal was located at the supersaturation where the compressibility did not reach a plateau over the course of the simulation. They found that the metastable region estimated in that way is much shallower as that calculated with the van der Waals equation of state and other mean field approaches. Here one should keep in mind that the compressibility increases when a droplet forms because one can easily compress a vapor if the system can evade by further condensation, i.e., the systems moves along the tie line, which is horizontal in the pressure-density diagram or its finite size counterpart. Therefore one can identify that spinodal by the onset of decomposition at given temperature consistent with

MD Simulations of Carbon Dioxide the picture of the so-called “effective spinodal” in finite size systems discussed above but not the thermodynamic mean field spinodal. Conclusions Thermodynamic stability limits are determined for pure carbon dioxide by molecular dynamics simulation based on a newly developed method1 using interfacial properties. The relation of the tangential pressure component through the vapor-liquid interface to the van der Waals loop in the nonstable region of the phase diagram is established for carbon dioxide. The correlation of the interface density profile with the tangens hyperbolicus and the error function gives similar agreement with the simulation data. The obtained thickness parameter is slightly bigger for the error-function profile because it inherently represents a wider interface. However these differences do not affect the estimation of the interface spinodal properties from the interface profiles. Since tangential pressure data cannot be measured experimentally, we also checked the applicability of eq 6 using a consistent set of simulation data for surface tension, interface thickness, and vapor pressure. The agreement between the spinodal estimated in this way to the spinodal obtained from the extreme values of the tangential pressure profile is quite good. This supports that eq 6 can be used to estimate liquid-vapor spinodal for pure liquids from experimentally accessible data. Furthermore, it follows that the shape factor s is unity for carbon dioxide as well as reported earlier for argon.1 The test of eq 6 with experimental data for helium 4 also supports this value.52 In conclusion the estimation of the liquid spinodal from experimental data is possible without adjustable parameter. The investigation of finite size effects on the estimation of the spinodal from interface properties reveals that the temperature-density spinodal data are not affected. A widening of the interface by finite size effects such as caused by capillary waves or any other influence affects the density profile and the tangential pressure profile simultaneously. Therefore, the extreme values of the tangential pressure are always at the same density values. Still there is an influence of the interface thickness on the spinodal pressure which requires further investigations. For small system sizes where capillary waves are not present, one obtains the mean field spinodal as obtained by an equation of state. Another important result of this investigation is that apparently very long simulation time periods are necessary to get reliable interface properties. Acknowledgment. The authors wish to explain their thanks to Dr. M. Moldover (NIST, USA) for his valuable comments and Dr. D. Paschek for friendly support with the program MOSCITO. Comments on the original manuscript by K. Binder are gratefully acknowledged. A. R. Imre was supported by Hungarian Research Fund (OTKA) under Contract No. K67930, the German Humboldt Foundation, and by the Bolyai Research Grant of the Hungarian Academy of Science. T. Kraska acknowledges support of the DFG by Grant KR1598/26-1. References and Notes (1) Imre, A. R.; Mayer, G.; Ha´zi, G.; Rozas, R.; Kraska, T. J. Chem. Phys. 2008, 128, 114708.

J. Phys. Chem. B, Vol. 113, No. 14, 2009 4697 (2) Binder, K. Phys. ReV. A 1984, 29, 341–349. (3) Volmer, M; Weber, A. Z. Phys. Chem. 1926, 119, 277–301. (4) Farkas, L. Z. Phys. Chem. 1927, 125, 236–248. (5) Becker, R.; Do¨ring, W. Ann. Phys. 1935, 24, 719–752. (6) Pinhasi, G. A.; Ullmann, A.; Dayan, A. Int. J. Heat Mass Transfer 2007, 50, 4780–4795. (7) Brunner, G. Supercritical Fluids as SolVents and Reaction Media; Elsevier: 2004. (8) Holloway, S. Philos. Trans. R. Soc. A 2007, 365, 1095–1107. (9) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021–12024. (10) Ciccotti, G.; Ferrario, M.; Ryckaert, J.-P. Mol. Phys. 1982, 47, 1253–1264. (11) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comp. Phys. 1977, 23, 327–341. (12) Chen, F.; Smith, P. E. J. Chem. Phys. 2007, 126, 221101. (13) Alejandre, J.; Tildesly, D. J.; Chapela, G. A. J. Chem. Phys. 1995, 102, 4574–4583. (14) Paschek, D.; Geiger, A. Moscito 4; University of Dortmund: 2002. (15) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (16) Leng, C. A.; Rowlinson, J. S.; Thompson, S. M. Proc. R. Soc. London, Ser. A 1976, 352, 1–23. (17) Weeks, J. D. J. Chem. Phys. 1977, 67, 3106–3121. (18) Mayer, G.; Hazi, G.; Pales, J.; Imre, A. R.; Fischer, B.; Kraska, T. Int. J. Mod. Phys. C 2004, 15, 1049–1060. (19) Buff, F. P.; Lovett, R. A.; Stillinger, F. H, Jr. Phys. ReV. Lett. 1965, 15, 621–623. (20) Wu, E. S.; Webb, W. W. Phys. ReV. A 1973, 8, 2065–2076. (21) Baeglehole, D. Phys. ReV. Lett. 1987, 58, 1434–1436. (22) Sengers, J. V.; van Leeuwen, J. M. J. Phys ReV A. 1989, 39, 6346– 6355. (23) Binder, K.; Mu¨ller, M. Int. J. Mod. Phys. C 2000, 11, 1093–1113. (24) Iring, J. H.; Kirkwood, J. G. J. Chem. Phys. 1950, 18, 817–829. (25) Bakker, G. Kapillarita¨t und Oberfla¨chenspannung, Handbuch der Experimentalphysik; Wien, W., Harms, F., Lenz, H., Eds.; Akad. Verlagsges: Leipzig, 1928. (26) Fuchs, K. Exner’s Rep. Phys. (Mu¨nchen) 1888, 24, 141–160. (27) van der Waals, J. D. Z. Phys. Chem. 1894, 42, 655–725. (28) Peng, D.-Y.; Robinson, D. P. Ind. Eng. Chem. Fundam. 1976, 15, 59–64. (29) Schmidt, J. W.; Moldover, M. R. J. Chem. Phys. 1993, 99, 582– 589. (30) Leonhard, K.; Kraska, T. J. Supercrit. Fluids 1999, 16, 1–10. (31) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon Press: 1989. (32) Orea, P.; Lo´pez-Lemus, J.; Alejandre, J. J. Phys. Chem. 2005, 123, 114702. (33) Gonza´lez-Melchor, M.; Bresme, F.; Alejandre, J. J. Chem. Phys. 2005, 122, 104710. (34) Smit, B. J. Phys. Chem. 1992, 96, 8639–8640. (35) MacDowell, L. G.; Virnau, P.; Mu¨ller, M.; Binder, K. J. Chem. Phys. 2004, 120, 5293–5308. (36) Shen, V. K.; Errington, J. R. J. Phys. Chem. B 2004, 108, 19595– 19606. (37) Kaski, K.; Binder, K.; Gunton, J. D. Phys. ReV. B 1984, 29, 3996– 4009. (38) Binder, K. Eur. Phys J. B 2008, 64, 307–314. (39) Gregory, V. P.; Schug, J. C. Mol. Phys. 1994, 82, 677–688. (40) Yamamoto, R.; Kitao, O.; Nakanishi, K. Mol. Phys. 1995, 84, 757– 768. (41) Kivelson, D.; Reiss, H. J. Phys. Chem. B 1999, 103, 8337–8343. (42) Hansen, J.; Verlet, L. Phys. ReV. 1969, 184, 151–161. (43) Corti, D. S.; Debenedetti, P. G. Chem. Eng. Sci. 1994, 49, 2717– 2734. (44) Nie, C.; Geng, J.; Marlow, W. H. J. Chem. Phys. 2008, 128, 234319. (45) Langer, J. S. Physica 1974, 73, 61–72. (46) Binder, K.; Stauffer, D. AdV. Phys. 1976, 25, 343–396. (47) Shneidman, V. A.; Jackson, K. A.; Beatty, K. M. Phys. ReV. B 1999, 59, 3579–3589. (48) Kraska, T. J. Chem. Phys. 2006, 124, 054507. (49) Yang, S. H.; Nosonovsky, M.; Zhang, H.; Chung, K.-H. Chem. Phys. Lett. 2008, 451, 88–92. (50) Werner, A.; Schmid, F.; Mu¨ller, M.; Binder, K. Phys. ReV. E 1999, 59, 728–738. (51) Heermann, D. W.; Wang, T.; Lee, W.-D. Int. J. Mod. Phys. C 2003, 14, 81–94. (52) Imre, A. R.; Kraska, T. Physica B 2008, 403, 3663–3666.

JP808789P

The Relation of Interface Properties and Bulk Phase Stability

Recommend Documents