Molecular Dynamics Simulation of Interaction between Supercritical

Jul 29, 2008 - The structural and dynamical properties of the supercritical CO2 fluid confined in the slit nanopores with the hydroxylated and silylat...
2 downloads 0 Views 2MB Size
J. Phys. Chem. C 2008, 112, 12815–12824

12815

Molecular Dynamics Simulation of Interaction between Supercritical CO2 Fluid and Modified Silica Surfaces Yan Qin,† Xiaoning Yang,*,† Yongfang Zhu,† and Jialun Ping‡ State Key Laboratory of Material-Orientated Chemical Engineering, College of Chemistry and Chemical Engineering, Nanjing UniVersity of Technology, Nanjing 210009, China, and Department of Physics, Nanjing Normal UniVersity, Nanjing 210097, China ReceiVed: December 20, 2007; ReVised Manuscript ReceiVed: June 11, 2008

The structural and dynamical properties of the supercritical CO2 fluid confined in the slit nanopores with the hydroxylated and silylated amorphous silica surfaces have been studied using molecular dynamics (MD) simulation. The amorphous bulk silica was obtained by a melt-quench MD simulation technique and the modified silica surfaces were artificially created by the attachment of hydrogen (-OH model) and trimethysilane (-Si(CH3)3 model) to the nonbridging oxygen atoms on the silica surfaces. The VdW interaction potential between the CO2 molecule and the hydroxylated silica surface was determined based on the ab initio quantum mechanics (QM) computation. The adsorption potential distributions of CO2 on the two modified silica surfaces were examined in order to evaluate the different surface interaction characteristics. The density profiles, the radial distribution functions, as well as the interfacial dynamics properties (self-diffusion coefficients and residence time) for the confined supercritical CO2 fluid have been simulated. It is demonstrated that the hydroxylated silica surface gives a stronger confining effect on the supercritical CO2 fluid as compared with the silylated surface. The remarkable impact on the supercritical CO2 fluid from the hydroxylated silica surface can be attributed to the H-bonding interaction between CO2 molecules and surface silanol groups. The analysis of the vibrational density of states of the confined supercritical CO2 fluid reveals the phenomena of the spectral shifts and the Fermi resonance in compaison with the bands in unconfined supercritical CO2. This spectrum behavior is associated with the enhanced interaction from the functional groups on silica surfaces. Introduction Supercritical CO2 (Sc-CO2) fluid has received considerable attention as an alternative to conventional solvents in the microelectronics industry.1–4 In particular, the Sc-CO2-based cleaning/drying technique is favorable to the highly confined geometry of advanced microelectronics materials/devices with nanoscale porosity.5 In the microelectronics field, amorphous silica is a key material owing to its special properties.6 A microscopic level understanding of the interaction between the Sc-CO2 fluid and the silica surface is very important in the development of the Sc-CO2 fluid applications. The study of the interaction between the CO2 fluid and the silica surface is also required in the chromatography/adsorption separation,7,8 and the heterogeneous catalyzed reaction9 in ScCO2 fluid. In these respects, porous silica materials with defined pore sizes usually serve as adsorbents, catalysts, and supports. Many applications of silica materials depend on their surface properties,10,11 which are largely determined by the functional groups on surfaces through various modification methods. In modified silica surfaces, the hydrophilic and hydrophobic surfaces are the main types for different purposes. It has been known that the silica surface reacts rapidly with moisture, leading to the formation of surface silanol groups, which possesses a hydrophilic feature. The concentration, distribution, and nature of silanols determine the related properties of the hydroxylated silica surface.12 The surface silylation have been * To whom correspondence should be addressed. E-mail: Yangxia@ njut.edu.cn. Tel: 86-025-83587185. † Nanjing University of Technology. ‡ Nanjing Normal University.

prepared by alkylsilane groups,1,13 which can be used to obtain the hydrophobic silica surfaces. Current experimental methods are usually difficult to clarify the above-mentioned microscopic level interaction mechanism,14–16 whereas molecular simulation can offer an effective means of investigating the microscopic structural and dynamical properties of confined fluids near solid surfaces. These interfacial properties reflect the fluid-solid interaction. Lee and Rossky17 studied the structural and dynamical properties of water on artificial hydrophilic silica surface. The structure behavior of organic liquids, such as acetone18 and methanol19 confined in cylindrical modified silica pores, has been investigated by molecular simulation. These simulation reports indicate that there exists an obvious H-bonding interaction between the surface hydroxyl group and these confined fluids. The adsorption of acetone/ water20 mixture in mesoprous silica modified by surface silylation was studied by using the Grand canonical Monte Carlo (GCMC) simulation. The result shows that the highly silylated mesoporous silica membranes may be effective on the mixture separation due to the hydrophobic effect. In addition, other new storage materials such as SiC nanotubes21 and modified carbon nanotubes22 have been reported to give an enhanced hydrogen adsorption due to the special modified surface feature. These cited research reports advise that surface structure modification can produce a unique interface interaction, which is the crucial influencing factor determining confined fluid properties. Some molecular simulation studies have been conducted on the confined behavior of CO2 molecules in the silica nanopores.23–25 However, relatively fewer studies focus on the modified silica surfaces, even though the interfacial behavior of CO2 fluid on the pristine silica surfaces has been reported.26,27 On amorphous

10.1021/jp711964e CCC: $40.75  2008 American Chemical Society Published on Web 07/29/2008

12816 J. Phys. Chem. C, Vol. 112, No. 33, 2008

Qin et al.

Figure 1. (a) Configuration of -OH surface; (b) configuration of -Si(CH3)3 surface. Atoms in the amorphous silica substrate are displayed in red (oxygen) and yellow (silicon). The purple represents the surface oxygen connected with hydrogen (white) or trimethysilane group. The trimethysilane groups are shown in green (silicon) and black (methyl group). (c) The distribution contour of adsorption energy on the -OH surface; (d) same as panel c but on the -Si(CH3)3 surface. Energy scale is shown on the right color bar.

silica surfaces, the tetrahedral network is distorted. Moreover, the addition of surface functional-groups enhances the surface heterogeneity. CO2 molecule possesses a higher quadrupole moment, thus, it could produce an enhanced interaction with the functional groups on silica surfaces,28,29 especially the hydroxyl group.30 Bakaev et al.31,32 obtained the adsorption isotherms of CO2 in silica surface by the GCMC method, showing that the CO2 adsorption is a sensitive probe of the atomic structures of amorphous silica surfaces. Since the confined behavior of Sc-CO2 fluid near the silica surfaces is still far from being well understood, it is very interesting to investigate the effect of different functional groups of silica surfaces on the interfacial structural and dynamical properties of CO2 molecules, namely, how much the Sc-CO2 fluid is significantly perturbed by the silica surfaces and how the surface functional groups interact with the CO2 molecule. This result will not only be helpful to identify the interaction mechanism between the modified silica surface and CO2 molecules but also be beneficial for controlling the adsorption of the solutes of interest from the Sc-CO2 fluid onto solid surfaces. In this work, MD simulations were carried out for the ScCO2 fluid confined between two parallel modified silica surfaces. The classical melt-quenching method33,34 was applied to obtain the amorphous bulk silica. We considered two kinds of modified surface models: hydroxylated and silylated silica surfaces. The influence of different modified silica surfaces on the structural and dynamical behavior of confined CO2 fluids was analyzed. For the silica materials, the confined behavior of fluid near the amorphous structure is more complicated than that near the ordered crystalline one.35–40 The modified crystalline β-cristobalite surfaces were also studied for comparison. In the following section, we describe the Models and Simulation

Details in this work, followed by a presentation of the Results and Discussions. The last section provides the Conclusions. Models and Simulation Details Surface Models. In the first step, a bulk amorphous silica solid (30.1 Å × 30.1 Å × 30.1 Å) with a density of 2.2 g/cm3 was obtained through the well-established classical melt-quench procedure.33,34 Then, the amorphous silica surface was relaxed by a 10 ps MD simulation34 with the bulk layer being kept immobilized. The detailed generation procedure for the silica substrate and surfaces was given in the Supporting Information. On the silica surface, we distinguished the nonbridging oxygen (NBO) atom that was bonded to less than two silicon atoms from the bridging oxygen (BO) atom in the bulk silica. The NBO density on the obtained silica surface is approximately 1.9/nm2, which is in good agreement with the report41 on vitreous silica surfaces. In the second step, based on the obtained amorphous silica, we artificially generated two modified surface models. The hydrogen atoms were directly connected to the NBOs on the amorphous surface, forming the so-called silanol group on the silica surface. The trimethysilane (-Si(CH3)3) groups were attached to the surface NBOs. A geometry optimization was performed to avoid an overlap between the surface groups and the surrounding atoms. The two modified surface models correspond to the amorphous silica surfaces saturated with the hydroxyl groups and the trimethysilane groups, which are denoted as the hydrophilic -OH (hydroxylated) surface model and the hydrophobic -Si(CH3)3 (silylated) surface model, respectively, in this work. Figure 1 shows the simulation configurations for the two modified surface models. The ordered silica surface was generated by fracturing the β-cristobalite structure on the [010] crystallographic face,35 since

Supercritical CO2 Fluid and Modified Silica Surfaces

J. Phys. Chem. C, Vol. 112, No. 33, 2008 12817

β-cristobalite is the crystalline phase of silica with the density and refractive index closest to the amorphous silica. The modified ordered surfaces were also created by connecting the hydrogen atoms and the trimethylsilane groups to the surface oxygen atoms. The total number of the trimethylsilane groups attached on the crystalline surface is half as many as the number of the hydroxyl groups due to the larger structure size for -Si(CH3)3. Figure S2 in Supporting Information shows the configuration illustrations of the modified crystalline surfaces. Intermolecular Potentials. In this simulation, the intermolecular interaction is described by the site-site Lennard-Jones potential and the coulomb contribution

Einter )

qiqj + 4εij[(σij/rij)12 - (σij/rij)6] rij

(1)

where σij and εij are the size and energy parameters, and qi is the charge of site i. For the CO2 molecule, a three-site model with a partial charge on each atom was used. The corresponding VdW parameters and the site charges were directly adopted from the EPM2 model.42,43 However, unlike the previous EPM2 model, the C-O bond was controlled to be flexible with the bond force constant kb in this work.44 The detailed CO2 model used in this work is given in the Supporting Information (Table S1). There are a large number of empirical force fields for SiO2. However, it is difficult to determine the reasonable interactive potential between CO2 molecules and SiO2 surfaces. For the -OH surface model, quantum mechanics (QM) computation was applied to obtain the mixed VdW parameters. Since the real silica surfaces are very complicated, a simple model cluster, which includes the hydroxyl groups, was designed to carry the most significant characteristics of the hydroxylated silica surface. Accordingly, the smaller cluster, (H3SiO)2Si(OH)2, including the SiO2 unit and silanol groups, was chosen as the model cluster structure (Supporting Information, Figure S3). This cluster structure has been successfully used to determine the force field parameters of the silica-water interface.45 The geometries and the atomic charges of the silica cluster model were optimized using the B3LYP functional at the 6-31G* basis set. The relevant geometric parameters are listed in Supporting Information, Table S2, which agree well with the structure characteristics of the model cluster reported in the previous references.45–47 In the determination of the mixed CO2-silica VdW interaction, we first obtained a series of interaction energies of various configurations between the optimized cluster and a CO2 molecule from the ab initio QM computations, and then these interaction energy data were fitted to the analytical potential function (eq 1) including the site-site Lennard-Jones+Coulomb interaction terms in order to fit the Lennard-Jones parameters between the silica cluster and CO2 molecule. The MP2/6-31G** level was used to obtain the interaction energy (see Supporting Information). Eighty-seven different interactive configurations were performed for the system. Four representative configurations are shown in Supporting Information, Figure S3. In the fitting of the VdW potential parameters for the silica cluster, only the silicon and oxygen atoms were considered and the parameter of the hydrogen atom in the hydroxyl group was directly adopted from the reference 45. The interaction energies from the 87 different configurations were fitted simultaneously to eq 1 by minimizing this following objective function

∑ 87

OBF )

(EC-C - EQM)2

(2)

i)1

where EC-C is the interaction energy between the cluster and CO2 molecule, and EQM is the corresponding energy from QM computation. The detailed fitting procedure was given in Supporting Information. The obtained mixed VdW parameters between CO2 and the hydroxylated silica surface are listed in Table 1, together with the initial parameter values. It should be emphasized that the fitted parameters are not unique. However, the choice of the initial values in the fitting procedure is based on the previous silica model. Hence, the obtained mixing VdW parameters of CO2-silica should be within the reasonable value range. In comparison with the starting parameters in the fitting, the values of σ have decreased slightly for the O-C(CO2) and O-O(CO2) interactions and the corresponding ε values have increased slightly. For the Si-O(CO2) and Si-C(CO2) interactions, the most significant changes are in the σ parameters. For the -Si(CH3)3 surface models, the interaction with CO2 molecules is only represented by the Lennard-Jones potential and only the oxygen atoms in the silica skeleton structure and the united CH3 groups on the surface are included in the interaction potential. The corresponding parameters are available directly from the reference 20 and they are listed in Table S3 of the Supporting Information. Simulation Details. In the simulation systems, two parallel silica surfaces obtained above were set to face each other in the Y direction. We selected the surface area of 30.1 × 30.1 Å2 in XZ plane for the amorphous surfaces, in which the thickness of the -OH surface is 16.8 Å in Y direction and 19.3 Å for the -Si(CH3)3 surface. The size of modified crystalline slab is 27.75 × 12.53 × 27.75 Å3. Except for the surface silanol groups, the whole silica slabs, including the -Si(CH3)3 groups, were kept rigid during the simulation process. For the hydroxylated surface, the silanol structure was flexible and governed by the harmonic potentials with the bond and angle force constants (see Supporting Information, Table S4). In this study, two separations (54, 20 Å) were investigated between the two parallel surfaces. The separation is defined as the distance between the outmost surface atoms/groups in the slit pores. The Sc-CO2 fluid phase, obtained from the bulk-phase MD simulations, was inserted into the confined slit spaces. Since the amorphous surfaces have irregular structure, it is not easy to determine precisely the slit pore volumes and the confined fluid densities. Thus, in this work, the available pore volumes for confined CO2 were estimated based on the distances between the average sites of NBO on the two parallel surfaces. The densities of confined CO2 fluid in the slit pores were approximately 0.650 and 0.225 g/cm3, respectively. For the crystalline surfaces, we only considered the case of large slit pore width. The site charge of the oxygen atom was determined from the optimized oxygen charge in the cluster structure and the charge of silicon atom was fixed with TABLE 1: CO2-Hydroxylated Silica Interaction Parameters type

fitting parameters σ (Å) ε (kcal/mol)

hydroxylated silica-CO2 Si-C(CO2) 2.607 Si-O(CO2) 2.745 O-C(CO2) 2.902 O-O(CO2) 3.043

0.172 0.291 0.114 0.194

initial parametersa σ (Å) ε (kcal/mol) 3.338 3.476 2.955 3.093

0.183 0.309 0.092 0.156

a The initial parameter values were obtained by the L-B mixing rule calculation between the CO2 model and reported hydroxylated silica model.45

12818 J. Phys. Chem. C, Vol. 112, No. 33, 2008

Qin et al.

Figure 2. Density distribution and free energy distribution of CO2 molecules in the slit pores between amorphous silica surfaces (54 Å). (a) F ) 0.65 g/cm3 for -OH surfaces; (b) F ) 0.65 g/cm3 for -Si(CH3)3 surfaces; (c) F ) 0.225 g/cm3 for -OH surfaces; (d) F ) 0.225 g/cm3 for -Si(CH3)3 surfaces. The dashed lines mark the average positions of hydrogen atoms and methyl groups.

qO ) -1/2qSi.48 The charge of hydrogen atom in the surface hydroxyl was adjusted to meet the neutral condition in the silica surfaces, with 0.414 e for the amorphous surfaces and 0.535 e for the ordered surfaces. The confined behavior of the CO2 fluid in the slit pores was studied using the classical MD simulations. All MD simulations were performed using the periodic boundary conditions in the NVT ensemble at 318.15 K. The Berendsen algorithm was employed to keep the desired temperature. The long-range coulomb interactions were treated using the Ewald summation algorithm, and a cutoff distance of 10 Å was applied to the VdW interactions for all cases. The MD equations of motion were integrated using the Beeman algorithm with a 1 fs time step. Simulation configuration sets were saved every 100 steps for data analysis. The total simulation duration for individual system was about 2.0 ns, in which the initial 1.0 ns simulation was used for the system equilibrium. Results and Discussions Adsorption Energy Distributions. In order to obtain a pictorial description of the interaction between the CO2 molecule and the modified amorphous silica surfaces, the adsorption energy distribution of a CO2 molecule over the surfaces was calculated by following the previous scheme,26,27 in which the minimumal fluid-solid potential energy was obtained for each grid point in the surfaces. The obtained results are shown in Figure 1c,d as an energy contour plot for the two modified silica surfaces. It is clear that the -OH surface possesses stronger interaction sites, where the mean potential energy is in the range of 5-8 kcal/mol with the lowest potential of 12 kcal/mol. The mean potential energy is compatible with the dilution adsorption heat (in range of 4.8-8.4 kcal/mol) of CO2 measured in actual silica solids.49,50 In contrast, the -Si(CH3)3 surface shows relatively weaker interaction energy with the mean potential being 3-5 kcal/mol. The comparative result indicates that the -OH surface could produce stronger attraction for CO2 molecules than the -Si(CH3)3 surface. In Figure 1, the potential distribution does not show a consistent pattern with the surface structures (Figure 1a,b). For example, in the -OH surface the

Figure 3. Density distribution of CO2 molecules in the slit pores between modified crystalline silica surfaces (54 Å). (a) -OH surfaces; (b) -Si(CH3)3 surfaces. The dashed lines denote the positions of the outmost hydrogen atoms and methyl groups.

lower potential energy sites do not exactly correspond to the locations of the surface hydroxyl groups. This may be due to the amorphous feature on the surface, where the various functional groups (Si-O-H and Si-O-Si) with the distorted tetrahedron structure jointly generate the enhanced interaction with the CO2 molecule. Structure Properties. The one-dimensional density profiles of the oxygen and carbon atoms of CO2 molecules in the slit pores of 54 Å with the modified amorphous silica surfaces are shown in Figure 2. The blue dashed lines roughly mark the average positions of the hydrogen atoms in the surface hydroxyl groups. Please note that for the -Si(CH3)3 surfaces the line indicates the average position of the CH3 groups. It is clearly observed that interfacial CO2 molecules appear evident overlap with the surface silanol layers, meaning that the CO2 molecules could, to great extent, access the amorphous -OH surfaces. However, on the Si(CH3)3-modified surfaces the interpenetration extent of CO2 molecules into the trimethylsilane layer is small.

Supercritical CO2 Fluid and Modified Silica Surfaces

J. Phys. Chem. C, Vol. 112, No. 33, 2008 12819

Figure 4. Two-dimensional density (in g/cm3) distribution contour of CO2 molecules on the modified silica surfaces with density F ) 0.65 g/cm3 for large distance separations (54 Å): (a) the hydroxylated amorphous surface; (b) the silylated amorphous surface. Density scale is shown on the right color bar.

Near the amorphous silica surfaces, the density profiles of CO2 display wide and smooth peaks with weak layering phenomenon. This observation has been reported in the previous simulations on the pristine silica surfaces26 and the amorphous MCM-41 pore surface.51 As shown in Figure 2a, the -OH surface induces a higher CO2 density peak and a closer peak position with respect to the surface as compared with the -Si(CH3)3 surface, near which a lower density peak with a farther position is observed. This difference is associated with the various interaction strengths of the two silica surfaces (see Figure 1). For the condition of low density (Figure 2c,d), an analogous behavior is observed for the two types of surfaces with the exception of obvious CO2 depletion in the central region of the pores. In the case of small slit pores of 20 Å, a similar result is obtained (Supporting Information, Figure S5). As a whole, the hydroxylated silica surface demonstrates a stronger CO2-philic feature, whereas the silylated surface has weaker attraction for CO2 molecules. Figure 3 shows the CO2 (at F ) 0.65 g/cm3) density profiles near the modified crystalline surfaces. The dashed lines also denote the positions of the outmost surface atoms/groups. For the crystalline counterparts, the -OH surface exhibits a similar strong attraction for CO2 molecules, whereas CO2 molecules are remarkably accumulated far away from the silylated surface, with about 3 Å depletion region from the Si(CH3)3 layer. Moreover, unlike the amorphous surfaces, distinct sharp density peaks are found near each crystalline surface (Figure 3). The layering behavior on the ordered silica surfaces may be as a result of the stronger regular interaction from the well-organized surface structures. The interfacial accumulation of CO2 molecules can be clearly represented by the two-dimensional density distribution in Figure 4, where the CO2 molecules in the interfacial region (27 Å < |Y| < 32 Å near the amorphous surfaces) were averaged along the Y direction as a function of the parallel cross-section positions. For the silanol surfaces (Figure 4a), the interfacial CO2 molecules show noticeable inhomogeneous distribution features, that is, the CO2 molecules penetrate into the interstices of the surface hydroxyl groups and locate near or around the hydroxyl groups with separated clusters. This configuration arrangement allows the CO2 molecules to interact from the side direction with the surface hydroxyl groups. On the other hand, relatively fewer CO2 molecules can accumulate on the silylated surface. Careful inspection of Figure 4b indicates that CO2 molecules can access the surface mainly through the special interfacial region where fewer functional groups (Si(CH3)3) appear. For the modified crystalline silica surfaces, the twodimension density distributions are given in Figure S6 (in the

Figure 5. The in-plane radial distribution functions of CO2 molecules (F ) 0.65 g/cm3) between the modified surfaces for large distance separation (54 Å). (a) Amorphous -OH surface; (b) amorphous -Si(CH3)3 surface; (c) crystalline -OH surface; (d) crystalline -Si(CH3)3 surface.

Supporting Information), which gives an additional verification of different interfacial interactions between the two modified silica surfaces and CO2 molecules. The free energy profiles of a CO2 molecule near the amorphous surfaces can be evaluated through a change in the free energy of a CO2 molecule from the bulk phase to the interface region. The free energy was calculated by integrating the average force of the CO2 molecule in the direction perpendicular to the surface from the center region to the surfaces over 1 ns MD run.52,53 Figure 2 and Supporting Information, Figure S5, illustrate the free energy profiles of a CO2 molecule within the slit pores. It is observed that the free energy and the density profiles show the reasonable “mirror symmetry” feature.26 The pattern of the free energy distribution suggests that the transfer of a CO2 molecule from the bulk region to the solid surface is thermodynamically favorable. The free energy near the -OH surface is lower than that on the -Si(CH3)3 surface, which is in accordance with the different surface affinity to CO2 molecules. The in-plane C-C radial distribution functions (RDF) of CO2 molecules in the confined slit pores have been given in Figure 5 for different interface regions. These curves were obtained by dispensing with the Y-coordinate in the distance calculation. On the basis of the confined CO2 locations, three slabs were defined (see Supporting Information, Table S5): adsorbed slab for CO2 molecules in the contacting layer; central slab in the

12820 J. Phys. Chem. C, Vol. 112, No. 33, 2008

Figure 6. (a) Intermolecular RDFs between the surface hydroxyl groups and the interfacial CO2 molecules. Red lines: for the crystalline surface. Black lines: for the amorphous surface. The inset shows the distribution of the angle O(CO2)...O(OH)-H(OH) for CO2 adsorbed on amorphous surfaces. (b) The typical adsorption configuration of CO2 molecules around silanol groups on silica surface. The capsule represents the CO2 molecule. The silanol groups are shown in purple (oxygen) and white (hydrogen).

middle region of the slit pores; intermediate slab for the region between the former two slabs. The positions and intensities in the RDFs of the intermediate and central slabs for the amorphous surfaces are almost identical to each other, and the RDF in the adsorbed layer shows somewhat higher first peak for the -OH surface. On the whole, the presence of the two modified amorphous silica surfaces does not cause apparent distortion in the fluid structure of the confined Sc-CO2 fluid. This result is somewhat different from the observation on the structured surfaces, such as the clay54 and the regular graphite55 surfaces. For comparison, the in-plane RDFs of CO2 molecules on the hydroxylated ordered silica surface have also been shown in Figure 5, where the RDF in the contact layer of the ordered silanol surface exhibits obvious ordered structure feature (Figure 5c). According to the above analysis, the hydrophilic hydroxylated surfaces produce stronger interaction for CO2 molecules than the hydrophobic silylated surfaces. This interaction mechanism can be examined by the RDFs of the interfacial CO2 molecules with respect to the surface hydroxyl groups. Figure 6a shows the RDFs between the surface hydroxyl groups and the oxygen atoms (CO2) in the adsorbed slabs for both the amorphous and the ordered hydroxylated surfaces. The first sharp peak position in the O(CO2)-O(hydroxyl) RDF is around 3.09 Å on the amorphous surface and 2.91Å on the crystalline surface. The corresponding O(CO2)-H(hydroxyl) RDFs also illustrate some characteristic peaks in the positions of 2.06 Å. However, the RDFs of the CO2 molecules with regard to the surface bridging oxygen (BO) atoms (in SiO2) do not show a noticeable

Qin et al. structural feature. The results indicate that the CO2 molecules in the contacting layers statistically present a preferentially shortrange order around the surface hydroxyl groups. The inset of Figure 6a displays the possibility distribution of the forming angle ∠O(CO2)...O(OH)-H(OH) on the amorphous -OH surface. There appears a bimodal distribution for the angle. The typical adsorption configurations of CO2 molecules around silanol groups are shown in Figure 6b. This geometrical arrangement (mode I) between the CO2 molecules and the surface silanol groups is typical of hydrogen-bonding formation56 and is consistent with the RDFs in Figure 6a. Mode II represents another contacting style between the adsorbed CO2 molecules and surface silanol groups. From the definition of H-bonding,56 the ratio of the silanol groups on which the CO2 molecules are adsorbed as H-bonding is approximately 18% for the amorphous surface. Other silanol groups have the non H-bonding interaction with CO2 molecules. However, based on the RDF of O(CO2)-O(hydroxyl) for the amorphous surface in Figure 6a, the ratio of the silanol groups, around which the CO2 molecules are located within the first peak position, is about 50%. This statistical analysis suggests that the surface silanol groups present an obviously localized contact with CO2 molecules, which also makes a contribution to the strong interaction between CO2 molecules and the hydroxylated surface. This possible H-bonding formation can be further verified by the additional QM optimization calculation between the SiO2 model cluster and a CO2 molecule. The optimization geometry of the model system (Supporting Information, Figure S7) captures the H-bonding formation information (for details, see Supporting Information). It is universally accepted that the silica surface with silanol groups can form H-bonds with water,57 methanol19 and other organic molecules.58,59 In addition, the experimental study60 has demonstrated the existence of Hbonding between CO2 and the silica surface modified by alkylamino groups. In this simulation, the obvious stronger adsorption affinity for CO2 molecules on the silanol surface should be to a large extent ascribed to this specific H-bonding interaction between CO2 molecules and the surface hydroxyl groups. Recently, Chaffee30 observed a clear tendency for CO2 molecules to form H-bonding with the hydroxyl groups on hexagonal mesoporous silica surface, which is in good agreement with our result. Dynamical Properties. To gain insight into the effect of the modified silica surfaces on the diffusion behavior of interfacial CO2 molecules, we calculated the CO2 self-diffusion coefficient using the Einstein relation

Dxz ) lim 〈∆x2(t) + ∆z2(t)〉/(4t) t f∞

Dy ) lim 〈∆y2(t)〉/(2t)

(3)

t f∞

where Dxz and Dy are the diffusion coefficients parallel and perpendicular to the surfaces, respectively. In principle, the Einstein relation may not applicable to the diffusion calculation near solid surfaces. In this study, we adopted the scheme of Choudhury and Pettitt,61 in which only the molecules that remain in the region of interest at the initial time were considered. This scheme takes into consideration the crossing movement of molecules through the defined slab boundary in relatively long time interval, which could more reasonably represent the short-time dynamics behavior in confined regions. We calculated the perpendicular (Y) and parallel components (XZ) of the mean square displacements (MSDs) of CO2, as a function of the time, for different regions near the modified silica surfaces. Figure S8 in Supporting

Supercritical CO2 Fluid and Modified Silica Surfaces

J. Phys. Chem. C, Vol. 112, No. 33, 2008 12821 Ni



1 R(t) ) θ (t )θ (t + t) Ns i)1 i 0 i 0

Figure 7. The CO2 (F ) 0.65 g/cm3) diffusion coefficients versus the distance from the surfaces for both the -OH (red) and -Si(CH3)3 (black) surfaces for larger distance separation (54 Å). Dxz and Dy are the diffusion coefficients parallel and perpendicular to the surfaces, respectively.

Information shows the MSDs of CO2 at various interfacial regions for the hydroxylated silica surface. The linear relationship between the MSDs and time (