Molecular simulation study of cavity-generated instabilities in the superheated Lennard-Jones liquid

Previous equilibrium-based density-functional theory (DFT) analyses of cavity formation in the pure component superheated Lennard-Jones (LJ) liquid [S. Punnathanam and D. S. Corti, J. Chem. Phys. 119, 10224 (2003); M. J. Uline and D. S. Corti, Phys. Rev. Lett. 99, 076102 (2007)] revealed that a thermodynamic limit of stability appears in which no liquidlike density profile can develop for cavity radii greater than some critical size (being a function of temperature and bulk density). The existence of these stability limits was also verified using isothermal-isobaric Monte Carlo (MC) simulations. To test the possible relevance of these limits of stability to a dynamically evolving system, one that may be important for homogeneous bubble nucleation, we perform isothermal-isobaric molecular dynamics (MD) simulations in which cavities of different sizes are placed within the superheated LJ liquid. When the impermeable boundary utilized to generate a cavity is removed, the MD simulations show that the cavity collapses and the overall density of the system remains liquidlike, i.e., the system is stable, when the initial cavity radius is below some certain value. On the other hand, when the initial radius is large enough, the cavity expands and the overall density of the system rapidly decreases toward vaporlike densities, i.e., the system is unstable. Unlike the DFT predictions, however, the transition between stability and instability is not infinitely sharp. The fraction of initial configurations that generate an instability (or a phase separation) increases from zero to unity as the initial cavity radius increases over a relatively narrow range of values, which spans the predicted stability limit obtained from equilibrium MC simulations. The simulation results presented here provide initial evidence that the equilibrium-based stability limits predicted in the previous DFT and MC simulation studies may play some role, yet to be fully determined, in the homogeneous nucleation and growth of embryos within metastable fluids.     

Molecular dynamics in the isothermal-isobaric ensemble: the requirement of a "shell" molecule. II. Simulation results

The results of a series of constant pressure and temperature molecular-dynamics (MD) simulation studies based on the rigorous shell particle formulation of the isothermal-isobaric (NpT) ensemble are presented. These MD simulations validate the newly proposed constant pressure equations of motion in which a "shell" particle is used to define uniquely the volume of the system [M. J. Uline and D. S. Corti, J. Chem. Phys. (to be published), preceding paper]. Ensemble averages obtained with the new MD NpT algorithm match the ensemble averages obtained using the previously derived shell particle Monte Carlo NpT method [D. S. Corti, Mol. Phys. 100, 1887 (2002)]. In addition, we also verify that the Hoover NpT MD algorithm [W. G. Hoover, Phys. Rev. A 31, 1695 (1985); 34, 2499 (1986)] generates the correct ensemble averages, though only when periodic boundary conditions are employed. The extension of the shell particle MD algorithm to multicomponent systems is also discussed, in which we show for equilibrium properties that the identity of the shell particle is completely arbitrary when periodic boundary conditions are applied. Self-diffusion coefficients determined with the shell particle equations of motion are also identical to those obtained in other ensembles. Finally, since the mass of the shell particle is known, the system itself, and not a piston of arbitrary mass, controls the time scales for internal pressure and volume fluctuations. We therefore consider the effects of the shell particle on the dynamics of the system. Overall, the shell particle MD algorithm is an effective simulation method for studying systems exposed to a constant external pressure and may provide an advantage over other existing constant pressure approaches when developing nonequilibrium MD methods.     

Molecular dynamics in the isothermal-isobaric ensemble: the requirement of a "shell" molecule. III. Discontinuous potentials

Based on the approach of Gruhn and Monson [Phys. Rev. E 63, 061106 (2001)], we present a new method for deriving the collisions dynamics for particles that interact via discontinuous potentials. By invoking the conservation of the extended Hamiltonian, we generate molecular dynamics (MD) algorithms for simulating the hard-sphere and square-well fluids within the isothermal-isobaric (NpT) ensemble. Consistent with the recent rigorous reformulation of the NpT ensemble partition function, the equations of motion impose a constant external pressure via the introduction of a shell particle of known mass [M. J. Uline and D. S. Corti, J. Chem. Phys. 123, 164101 (2005); 123, 164102 (2005)], which serves to define uniquely the volume of the system. The particles are also connected to a temperature reservoir through the use of a chain of Nose-Hoover thermostats, the properties of which are not affected by a hard-sphere or square-well collision. By using the Liouville operator formalism and the Trotter expansion theorem to integrate the equations of motion, the update of the thermostat variables can be decoupled from the update of the positions of the particles and the momentum changes upon a collision. Hence, once the appropriate collision dynamics for the isobaric-isenthalpic (NpH) equations of motion is known, the adaptation of the algorithm to the NpT ensemble is straightforward. Results of MD simulations for the pure component square-well fluid are presented and serve to validate our algorithm. Finally, since the mass of the shell particle is known, the system itself, and not a piston of arbitrary mass, controls the time scales for internal pressure and volume fluctuations. We therefore consider the influence of the shell particle algorithm on the dynamics of the square-well fluid.     

On the equilibrium structures of the complexes H2C3H+ · Ar and c-C3H3(+) · Ar: results of explicitly correlated coupled cluster calculations

Explicitly correlated coupled cluster theory at the CCSD(T)-F12x (x = a, b) level [T. B. Adler et al., J. Chem. Phys. 127, 221106 (2007)] has been employed in a study of the potential energy surfaces for the complexes H(2)C(3)H(+) · Ar and c-C(3)H(3)(+) · Ar. For the former complex, a pronounced minimum with C(s) symmetry was found (D(e) ≈ 780 cm(-1)), well below the local "H-bound" minimum with C(2v) symmetry (D(e) ≈ 585 cm(-1)). The absorption at 3238 cm(-1) found in the recent infrared photodissociation spectra [A. M. Ricks et al., J. Chem. Phys. 132, 051101 (2010)] is, thus, interpreted as an essentially free acetylenic CH stretching vibration of the propargyl cation. A global minimum of C(s) symmetry was also obtained for c-C(3)H(3)(+) (D(e) ≈ 580 cm(-1)), but the energy difference with respect to the local C(2v) minimum is only 54 cm(-1).     

The free energy of the metastable supersaturated vapor via restricted ensemble simulations. II. Effects of constraints and comparison with molecular dynamics simulations

Extensive restricted canonical ensemble Monte Carlo simulations [D. S. Corti and P. Debenedetti, Chem. Eng. Sci. 49, 2717 (1994)] were performed. Pressure, excess chemical potential, and excess free energy with respect to ideal gas data were obtained at different densities of the supersaturated Lennard-Jones (LJ) vapor at reduced temperatures from 0.7 to 1.0. Among different constraints imposed on the system studied, the one with the local minimum of the excess free energy was taken to be the approximated equilibrium state of the metastable LJ vapor. Also, a comparison of our results with molecular dynamic simulations [A. Linhart et al., J. Chem. Phys. 122, 144506 (2005)] was made.     

Molecular dynamics in the isothermal-isobaric ensemble: the requirement of a "shell" molecule. I. Theory and phase-space analysis

Current constant pressure molecular-dynamics (MD) algorithms are not consistent with the recent reformulation of the isothermal-isobaric (NpT) ensemble. The NpT ensemble partition function requires the use of a "shell" molecule to identify uniquely the volume of the system, thereby avoiding the redundant counting of configurations [e.g., G. J. M. Koper and H. Reiss, J. Phys. Chem. 100, 422 (1996); D. S. Corti, Phys. Rev. E, 64, 016128 (2001)]. So far, only the NpT Monte Carlo method has been updated to allow the system volume to be defined by a shell particle [D. S. Corti, Mol. Phys. 100, 1887 (2002)]. A shell particle has yet to be incorporated into MD simulations. The proper modification of the NpT MD algorithm is therefore the subject of this paper. Unlike Andersen's method [H. C. Andersen, J. Chem. Phys. 72, 2384 (1980)] where a piston of unknown mass serves to control the response time of volume fluctuations, the newly proposed equations of motion impose a constant external pressure via the introduction of a shell particle of known mass. Hence, the system itself sets the time scales for pressure and volume fluctuations. The new algorithm is subject to a number of fundamentally rigorous tests to ensure that the equations of motion sample phase space correctly. We also show that the Hoover NpT algorithm [W. G. Hoover, Phys. Rev. A. 31, 1695 (1985); 34, 2499 (1986)] does sample phase correctly, but only when periodic boundary conditions are employed.     

Calculation of solid-liquid interfacial free energy: a classical nucleation theory based approach

We present a simple approach to calculate the solid-liquid interfacial free energy. This new method is based on the classical nucleation theory. Using the molecular dynamics simulation, we employ spherical crystal nuclei embedded in the supercooled liquids to create an ideal model of a homogeneous nucleation. The interfacial free energy is extracted by fitting the relation between the critical nucleus size and the reciprocal of the critical undercooling temperature. The orientationally averaged interfacial free energy is found to be 0.302+/-0.002 (in standard LJ unit). The temperature dependence of the interfacial free energy is also obtained in this work. We find that the interfacial free energy increases slightly with increasing temperature. The positive temperature coefficient of the interfacial free energy is in qualitative agreement with Spaepen's analysis [Solid State Phys. 47, FS181 (1994)] and Turnbull's empirical estimation [J. Appl. Phys. 21, 1022 (1950)].     

Bubble nucleation in polymeric liquids. II. theoretical considerations

A theory based on classical nucleation theory is developed for bubble nucleation in polymer solutions. The theory requires information on solubility, diffusivity, concentration, surface tension, temperature, and degree of supersaturation. The effects of supersaturation and of the presence of large molecules in a liquid mixture on the free energy of bubble formation are included in the theoretical development. A semiempirical equation for the determination of bubble nucleation rate is developed, with the aid of experimental results reported in part I of this series. Using the experimental data, computer simulations of bubble nucleation in polymer solutions are performed. The consumptions of the volatile component in a liquid mixture, due to bubble nucleation and subsequent growth, and the variation of bubble nucleation rate during the expansion process are included in the simulation of the bubble nucleation process.     

Monte Carlo simulation study of droplet nucleation

A new rigorous Monte Carlo simulation approach is employed to study nucleation barriers for droplets in Lennard-Jones fluid. Using the gauge cell method we generate the excess isotherm of critical clusters in the size range from two to six molecular diameters. The ghost field method is employed to compute the cluster free energy and the nucleation barrier with desired precision of (1-2)kT. Based on quantitative results obtained by Monte Carlo simulations, we access the limits of applicability of the capillarity approximation of the classical nucleation theory and the Tolman equation. We show that the capillarity approximation corrected for vapor nonideality and liquid compressibility provides a reasonable assessment for the size of critical clusters in Lennard-Jones fluid; however, its accuracy is not sufficient to predict the nucleation barriers for making practical estimates of the rate of nucleation. The established dependence of the droplet surface tension on the droplet size cannot be approximated by the Tolman equation for small droplets of radius less than four molecular diameters. We confirm the conclusion of ten Wolde and Frenkel [J. Chem. Phys. 109, 9901 (1998)] that integration of the normal component of the Irving-Kirkwood pressure tensor severely underestimates the nucleation barriers for small clusters.     

Microdroplet dissolution into a second-phase solvent using a micropipet technique: test of the Epstein-Plesset model for an aniline-water system

The Epstein-Plesset model was originally derived for the dissolution of a single gas bubble in an infinite aqueous solution (Epstein, P. S.; Plesset, M. S. J. Chem. Phys. 1950, 18, 1505-1509). The micropipet manipulation technique was previously shown to test this theory on air microbubbles and air-filled lipid-coated microparticles accurately and appropriately (Duncan, P. B.; Needham, D. Langmuir 2004, 20, 2567-2578). This same theory is now tested to model liquid microdroplet dissolution in a well-defined solution environment. As presented previously for the gas-bubble system, holding a single microparticle at the end of a micropipet was not shown to affect the dissolution profile and allowed isotropic diffusion significantly, a necessary condition for the validation of the theory. Here, an aniline-water system with an initial droplet diameter of 50 microm was used as a model liquid-liquid system. A microdroplet of aniline in an aqueous solution presatureated with aniline at distinct levels was tested, as was the reverse system of a water droplet in an aniline solution. The dissolution lifetime was shown to increase with increasing medium saturation fraction according to the Epstein-Plesset time-dependent theory (including the time required to establish the stationary layer) neglecting interfacial tension. The droplet lifetime can be increased by an order of magnitude (from about 10 to 100 s) by increasing the saturation fraction from 0 to 0.9 and by another order of magnitude by increasing from 0.9 to 0.99. The technique proved to be an accurate and appropriate method to test the dissolution of single liquid microdroplets in a second liquid solution and establishes a systematic experimental and theoretical approach to the investigation of the formation of polymer and other microparticles.     

Vapor-to-droplet transition in a Lennard-Jones fluid: simulation study of nucleation barriers using the ghost field method

We report a comprehensive Monte Carlo (MC) simulation study of the vapor-to-droplet transition in Lennard-Jones fluid confined to a spherical container with repulsive walls, which is a case study system to investigate homogeneous nucleation. The focus is made on the application of a modified version of the ghost field method (Vishnyakov, A.; Neimark, A. V. J. Chem. Phys. 2003, 119, 9755) to calculate the nucleation barrier. This method allows one to build up a continuous trajectory of equilibrium states stabilized by the ghost field potential, which connects a reference droplet with a reference vapor state. Two computation schemes are employed for free energy calculations, direct thermodynamic integration along the constructed trajectory and umbrella sampling. The nucleation barriers and the size dependence of the surface tension are reported for droplets containing from 260 to 2000 molecules. The MC simulation study is complemented by a review of the simulation methods applied to computing the nucleation barriers and a detailed analysis of the vapor-to-droplet transition by means of the classical nucleation theory.     

Dynamics of ice nucleation on water repellent surfaces

Prevention of ice accretion and adhesion on surfaces is relevant to many applications, leading to improved operation safety, increased energy efficiency, and cost reduction. Development of passive nonicing coatings is highly desirable, since current antiicing strategies are energy and cost intensive. Superhydrophobicity has been proposed as a lead passive nonicing strategy, yet the exact mechanism of delayed icing on these surfaces is not clearly understood. In this work, we present an in-depth analysis of ice formation dynamics upon water droplet impact on surfaces with different wettabilities. We experimentally demonstrate that ice nucleation under low-humidity conditions can be delayed through control of surface chemistry and texture. Combining infrared (IR) thermometry and high-speed photography, we observe that the reduction of water-surface contact area on superhydrophobic surfaces plays a dual role in delaying nucleation: first by reducing heat transfer and second by reducing the probability of heterogeneous nucleation at the water-substrate interface. This work also includes an analysis (based on classical nucleation theory) to estimate various homogeneous and heterogeneous nucleation rates in icing situations. The key finding is that ice nucleation delay on superhydrophobic surfaces is more prominent at moderate degrees of supercooling, while closer to the homogeneous nucleation temperature, bulk and air-water interface nucleation effects become equally important. The study presented here offers a comprehensive perspective on the efficacy of textured surfaces for nonicing applications.     

Curvature dependence of surface free energy of liquid drops and bubbles: A simulation study

We study the excess free energy due to phase coexistence of fluids by Monte Carlo simulations using successive umbrella sampling in finite L×L×L boxes with periodic boundary conditions. Both the vapor-liquid phase coexistence of a simple Lennard-Jones fluid and the coexistence between A-rich and B-rich phases of a symmetric binary (AB) Lennard-Jones mixture are studied, varying the density ρ in the simple fluid or the relative concentration x(A) of A in the binary mixture, respectively. The character of phase coexistence changes from a spherical droplet (or bubble) of the minority phase (near the coexistence curve) to a cylindrical droplet (or bubble) and finally (in the center of the miscibility gap) to a slablike configuration of two parallel flat interfaces. Extending the analysis of Schrader et al., [Phys. Rev. E 79, 061104 (2009)], we extract the surface free energy γ(R) of both spherical and cylindrical droplets and bubbles in the vapor-liquid case and present evidence that for R→∞ the leading order (Tolman) correction for droplets has sign opposite to the case of bubbles, consistent with the Tolman length being independent on the sign of curvature. For the symmetric binary mixture, the expected nonexistence of the Tolman length is confirmed. In all cases and for a range of radii R relevant for nucleation theory, γ(R) deviates strongly from γ(∞) which can be accounted for by a term of order γ(∞)/γ(R)-1∝R(-2). Our results for the simple Lennard-Jones fluid are also compared to results from density functional theory, and we find qualitative agreement in the behavior of γ(R) as well as in the sign and magnitude of the Tolman length.     

Thermodynamics and kinetics of bubble nucleation: Simulation methodology

The simulation of homogeneous liquid to vapor nucleation is investigated using three rare-event algorithms, boxed molecular dynamics, hybrid umbrella sampling Monte Carlo, and forward flux sampling. Using novel implementations of these methods for efficient use in the isothermal-isobaric ensemble, the free energy barrier to nucleation and the kinetic rate are obtained for a Lennard-Jones fluid at stretched and at superheated conditions. From the free energy surface mapped as a function of two order parameters, the global density and largest bubble volume, we find that the free energy barrier height is larger when projected over bubble volume. Using a regression analysis of forward flux sampling results, we show that bubble volume is a more ideal reaction coordinate than global density to quantify the progression of the metastable liquid toward the stable vapor phase and the intervening free energy barrier. Contrary to the assumptions of theoretical approaches, we find that the bubble takes on cohesive non-spherical shapes with irregular and (sometimes highly) undulating surfaces. Overall, the resulting free energy barriers and rates agree well between the methods, providing a set of complementary algorithms useful for studies of different types of nucleation events.     

The free energy of the metastable supersaturated vapor via restricted ensemble simulations

Pressure, excess chemical potential, and excess free energy, with respect to ideal gas data at different densities of the supersaturated Lennard-Jones particle vapor at the reduced temperature 0.7 are obtained by the restricted canonical ensemble Monte Carlo simulation method [D. S. Corti and P. Debenedetti, Chem. Eng. Sci. 49, 2717 (1994)]. The excess free energy values depend upon the constraints imposed on the system with local minima exhibited for densities below the spinodal density and monotonic variation for densities larger than the spinodal density. The results are compared with a molecular dynamics simulation study [A. Linharton et al., J. Chem. Phys. 122, 144506 (2005)] on the same system. The current study verifies the conclusion drawn by the simulation work that clustering of Lennard-Jones atoms exists even in the vicinity of spinodal. Our method gives an alternative to molecular dynamic simulations for the determination of equilibrium properties of a metastable fluid, especially close to the spinodal, and does not require a very large system to carry out the simulation.     

Relationships between the single particle barrier hopping theory and thermodynamic, disordered media, elastic, and jamming models of glassy systems

The predictions of the ultralocal limit of the activated hopping theory of highly viscous simple fluids and colloidal suspensions [K. S. Schweizer and G. Yatsenko, J. Chem. Phys. 127, 164505 (2007), preceding paper] for the relaxation time and effective activation barrier are compared with those of diverse alternative theoretical approaches and computer simulation. A nonlinear connection between the barrier height and excess pressure as empirically suggested by simulations of polydisperse repulsive force fluids is identified. In the dense normal and weakly dynamical precursor regime, where entropic barriers of hard spheres are nonexistent or of order the thermal energy, agreement with an excess entropy ansatz is found. In the random close packing or jamming limit, the barrier hopping theory predicts an essential singularity stronger than the free volume model, which is in agreement with the simplest entropic droplet nucleation and replica field theoretic approaches. Upon further technical simplification of the theory, close connections with renormalization group and nonperturbative memory function based studies of activated transport of a Brownian particle in a disordered medium can been identified. Several analytic arguments suggest a qualitative consistency between the barrier hopping theory and solid-state elastic models based on the high frequency shear modulus and a molecular-sized apparent activation volume. Implications of the analysis for the often high degeneracy of conflicting explanations of glassy dynamics are discussed.     

Ab initio molecular dynamics studies of tetrasulfur. Dynamics coupling the C2v open and D2h closed forms of S4

Ab initio molecular dynamics (AIMD) simulations were performed on the closed D(2h) and open C(2v) isomers of tetrasulfur. After a careful calibration of the electronic structure method, the calculations were done using the BPW91/aug-cc-pVTZ method. This combination of method/basis set adequately reproduces the relative benchmark CCSD(T) energy difference [Matus, M.; Dixon, D.; Peterson, K. A.; Harkless, J. A. W.; Francisco, J. S. J. Chem. Phys. 2007, 127, 174305] between these two isomers and, crucially, the fact that the D(2h) structure is a transition state linking two equivalent (mirror images) C(2v) isomers. The trajectories show that the symmetric open C(2v) isomers interconvert when passing through the D(2h) closed transition state structure and that, unlike tetraoxygen, no three-dimensional structures arise. The dynamic vibrational analysis yields peaks in good agreement with the static CCSD(T) harmonic frequencies and explains higher peaks as overtones, thus showing that unlike previous AIMD DFT-based approaches, carefully calibrated exchange-correlation functionals can produce reliable molecular dynamics results for complex PESs as the one corresponding to the lowest singlet of S(4).     

Ab initio treatment of the chemical reaction precursor complex Cl(2P)-HF. 2. Bound states and infrared spectrum

Bound energy levels and properties of the Cl(2P)-HF complex were obtained from full three-dimensional (3D) calculations, with the use of the ab initio computed diabatic potential surfaces from the preceding paper and the inclusion of spin-orbit coupling. For a better understanding of the dynamics of this complex we also computed a 2D model in which the HF bond length r was frozen at the vibrationally averaged values r0 and r1 and a 2 + 1D model in which the 3D potentials were averaged over the v(HF) = 0 and v(HF) = 1 vibrational wave functions of free HF. Also 1D calculations were made in which both r and the Cl-HF distance R were frozen. The complex is found to have the linear hydrogen bonded Cl-HF structure, with ground-state quantum numbers J = 3/2 for the overall angular momentum and /omega/ = 3/2 for its projection on the intermolecular axis R. The binding energy is D0 = 432.25 cm(-1) for v(HF) = 0 and D0 = 497.21 cm(-1) for v(HF) = 1. Bending modes with /omega/ = 1/2 and /omega/ = 5/2 are split by the Renner-Teller effect, since the electronic ground state is a degenerate 2pi state. A series of intermolecular (R) stretch modes was identified. Rotational constants and e-f parity splittings were extracted from the levels computed for J = 1/2 to 7/2. The computed red shift of the HF stretch frequency of 64.96 cm(-1) and the 35Cl-37Cl isotope shift of 0.033 cm(-1) are in good agreement with the values of 68.77 and 0.035 cm(-1) obtained from the recent experiment of Merritt et al. (Phys. Chem. Chem. Phys. 2005, 7, 67), after correction for the effect of the He nanodroplet matrix in which they were measured.     

Ab initio potential energy surfaces, bound states, and electronic spectrum of the Ar-SH complex

New ab initio potential energy surfaces for the (2)Pi ground electronic state of the Ar-SH complex are presented, calculated at the RCCSD(T)/aug-cc-pV5Z level. Weakly bound rotation-vibration levels are calculated using coupled-channel methods that properly account for the coupling between the two electronic states. The resulting wave functions are analyzed and a new adiabatic approximation including spin-orbit coupling is proposed. The ground-state wave functions are combined with those obtained for the excited (2)Sigma(+) state [D. M. Hirst, R. J. Doyle, and S. R. Mackenzie, Phys. Chem. Chem. Phys. 6, 5463 (2004)] to produce transition dipole moments. Modeling the transition intensities as a combination of these dipole moments and calculated lifetime values [A. B. McCoy, J. Chem. Phys. 109, 170 (1998)] leads to a good representation of the experimental fluorescence excitation spectrum [M.-C. Yang, A. P. Salzberg, B.-C. Chang, C. C. Carter, and T. A. Miller, J. Chem. Phys. 98, 4301 (1993)].