1764
J. Phys. Chem. C 2008, 112, 1764-1771
ARTICLES Molecular Dynamics Simulation Study of Growth Regimes during Polycondensation of Silicic Acid: from Silica Nanoparticles to Porous Gels Sudin Bhattacharya† Department of Mechanical Engineering, UniVersity of Michigan, Ann Arbor, Michigan 48109-2158
John Kieffer* Department of Materials Science and Engineering, UniVersity of Michigan, Ann Arbor, Michigan 48109-2136 ReceiVed: May 16, 2007; In Final Form: September 10, 2007
Molecular dynamics simulations based on a reactive force field with charge transfer were used to model the sol-gel synthesis of nanoporous silica gels in an aqueous environment. Three distinct growth regimes emerge, depending on the solvent-to-silica ratio: compact nanoparticles, percolated silica networks, and branched clusters. These growth regimes can be identified on the basis of distinctive structural features. In the case of compact particles, the radial distribution functions exhibit a broad maximum that coincides with the radius of gyration of the aggregates, whereas in continuous networks the radial distribution function increases steadily beyond the near-range structural features. Furthermore, these growth regimes can be distinguished on the basis of the concentrations of structural defects, such as dangling bonds and residual OH groups. The growth kinetics of individual regimes are characterized by different relative contributions of atomic diffusion to the overall aggregation rate. Finally, the resulting gel structures possess different mechanical stability, as can be assessed by quantifying the extents of structural collapse during simulated supercritical drying.
I. Introduction process1
The sol-gel allows for the preparation of ceramics and glasses with a wide range of structural and physical properties by controlling the chemical environment and processing parameters.2 There is thus considerable interest in understanding the mechanism of the sol-gel transformation and the features of the resulting structures, as reflected in a large number of studies in the literature, both experimental3-18 and computational.19-42 Sol-gel processing has been the method of choice for producing porous silica gels by the hydrolysis of silicon alkoxide precursors like TMOS [tetramethoxysilane, Si(OCH3)4] and TEOS [tetraethoxysilane, Si(OCH2CH3)4], followed by alcohol- and water-producing condensation.1 We focus in this article on simulating the water-producing condensation of aqueous silicic acid [Si(OH)4], that is, fully hydrolyzed precursor molecules, to yield polymerized silica gel structures as per the reaction:
≡Si-OH + HO-Si≡ f ≡Si-O-Si≡ + H2O
(1)
To simulate this process, we have used the molecular dynamics (MD) simulation technique based on a reactive threebody interatomic potential that allows for charge transfer upon the rupture or formation of chemical bonds.42-46 This approach makes it possible to simulate the elementary reaction processes * To whom correspondence should be addressed. E-mail: kieffer@ umich.edu, Phone: 734-763-2595, Fax: 734-763-4788. † Current address: Division of Computational Biology, The Hamner Institutes for Health Sciences, Research Triangle Park, NC 27709.
for sufficiently large numbers of atoms, with the objective of generating increasingly realistic gel structures that we can characterize in great detail. Previous MD simulations of sol-gel polymerization in silica systems include both reactive and nonreactive models. Pereira et al.37 have developed nonreactive potentials for liquids with varying compositions encountered in the sol-gel process water, methanol, ethanol, TMOS, and TEOS - at high temperature and pressure as well as at ambient conditions, based on prior ab initio studies.35,36 They have also applied these nonreactive potentials to indirectly model sol-gel reactions by simulating the constituents at three different stages of the process - before hydrolysis, after hydrolysis, and after condensation.38 The first working reactive MD simulations of sol-gel condensation were carried out by Feuston and Garofalini,21 using a modified Born-Mayer-Huggins (BMH) potential, with additional terms from the Rahman-Stillinger-Lemberg potential47 to describe hydrogen interactions. In a system consisting of 27 silicic acid monomers, the onset of polymerization was observed at an elevated temperature of 1500 K. This model was later extended by Garofalini and Martin to a system of 216 silicic acid monomers and longer run-times (120 ps) at temperatures of 525-650 K and relatively high densities (1.4-1.6 g/cm3).22 They observed the formation of chains in the early stages of polymerization, followed by ring formation - as confirmed by NMR studies and semiempirical quantum mechanical calculations. The time evolution of Qn silicon species (where the subscript n indicates the number of bridging oxygen atoms attached to a given silicon atom) was also found to be consistent
10.1021/jp073808f CCC: $40.75 © 2008 American Chemical Society Published on Web 01/19/2008
Polycondensation of Silicic Acid
J. Phys. Chem. C, Vol. 112, No. 6, 2008 1765
with NMR data. In a separate study based on the same model, Martin and Garofalini investigated the molecular mechanisms involved in the reaction and the effect of hydrogen.23 Rao and Gelb have recently carried out a large-scale reactive MD simulation41 of silicic acid polymerization in aqueous solvent based on the Feuston-Garofalini potential. Their model includes up to 729 Si(OH)4 monomers simulated for between 1.6 and 12.5 ns at high temperatures (1500-2500 K) and explores a broad range of water-to-silicon ratios and silicic acid concentrations. They observe cluster-cluster aggregation events leading to the formation of sol particles several nanometers in size, spanning up to 97% of the simulation box. They find that water promotes the formation of cyclic clusters and acts as a catalyst for the condensation reaction. The estimated activation energy of condensation (13-15 kcal/mol) in excess water is also found to be in good agreement with experimental data. In our previous work,42 we have successfully modeled solgel condensation with our reactive potential incorporating charge transfer and studied in detail the fractal dimension df of the resulting structures from geometric correlations. We also investigated the effect of supercritical drying on the fractal nature of condensed gels and compared their df values to those of ruptured silica glasses of comparable density. In the present study, we incorporate various amounts of additional water molecules into the initial reactive environment of our models (the only water molecules in our previous study were those that formed as a result of silicic acid condensation) to simulate more realistic chemical conditions and explore a wider range of system densities. We discuss the interrelation between the degree of reaction, the water-to-silicon ratio, ω, and the mobility of silicon atoms. We also examine the effect of density and ω on the fractal dimension of the condensed gels and the structures produced after supercritical drying. II. Computational Details To account for the mixed covalent-ionic nature of the Si-O bond in silica and its resultant charge localization and directionality, the potential model used in our simulations includes a Coulomb term incorporating charge transfer, a BMH repulsive term, a three-body term that accounts for covalent interactions, and a term to account for long-range dispersive interactions between nonbonded atoms. For a single atom, the potential energy is given by N
φi ) qi
qj
∑ j ) 1 4π r
0 ij NC - 1
NC
+
(
zj
i
j
∑ Aij 1 + n + n j)1
NC
∑ ∑
j)1 k)j+1
III. Results and Discussion
)
zi
tively, the BMH repulsive, the covalent attractive, and the van der Waals attractive components of the atomic interactions. The long-range effects of the Coulomb interactions, resulting from their dependence on the inverse of rij, are evaluated using the Ewald summation method.48,49 O-H interactions are also modeled using eq 2, that is, covalent contributions are described by the third term and hydrogen bonding by the fourth, appropriately parametrized. Further details about the individual terms and values of the specific parameters used for all of the pairs and three-body interactions involving silicon, oxygen, and hydrogen atoms have been provided in our previous publication.42 This potential was originally developed for the simulation of systems that exhibit mixed covalent-ionic bonding and has been successfully used to study structural transitions and anomalous thermomechanical properties in silica44,45,50-53 and boron oxide.54 We simulated the condensation of silicic acid molecules both in dry (i.e., no water initially present in the system) and aqueous environments. The dry configurations contained 512 silicic acid molecules (a total of 4608 atoms) in a cubic simulation box, with periodic boundary conditions applied. The system was relaxed at a temperature of 300 K, then heated to a high temperature (1000 K), and quenched to ensure that the molecules were dispersed randomly in the simulation box. By adjusting the size of the simulation box, we generated 9 different configurations at densities ranging from 0.2 to 1.8 g/cm3. For simulations in an aqueous environment, an ensemble of 216 silicic acid molecules and variable quantities of water molecules were used. The ratio of water to silicon atoms, ω, ranged from 1 to 5, and the resulting systems contained between 2592 and 5184 atoms. Corresponding to each water-to-silicon ratio, configurations at 9 different densities (0.2 to 1.8 g/cm3) were generated as above, resulting in a total of 54 systems. Condensation reactions were allowed to proceed at a temperature of 700 K to ensure reasonable simulation times (the reaction rate shows a steady increase between 300 and 700 K). The systems were typically allowed to react for between 1 and 3.5 ns. The equations of motion were integrated using a time step of ∆t ) 0.001 ps. Water molecules formed as a result of the condensation reaction were initially allowed to stay in the system - thus resembling the actual production of aerogels in autoclaves under high temperature and pressure conditions.
e(σi+σj-rij)Fij +
(φij + φik)e-γijk(θh -θijk) 2
N
Cij
j)1
rij6
∑
(2)
In the first (Coulomb) term, qi denotes the charge on the ith particle, 0 is the permittivity of free space, and rij is the distance between the ith and jth particles. The charge qi is evaluated as qi NC 0 ) q0i - 2 ∑j ) 1 δijζij, where qi is the charge of the isolated atom, NC is the coordination number of the atom when bonded, δij is the amount of transferred charge, and ζij ) (1 + e(rij-a)b)-1 is a function that modulates the transfer of charge with interatomic distance, a and b being empirical parameters. This provision of charge transfer allows for the formation and rupture of bonds in the reaction. Electroneutrality is maintained by having δji ) -δij. Covalent bonding is described by φij ) -Cij (κij/ηij) ζije(λij-rij)ηij, which is also modulated by the chargetransfer term. The remaining terms in eq 2 represent, respec-
We begin with a discussion of the observed condensation rates. The condensation reaction produces a polymerized silica network and water molecules: the degree of reaction can be measured as the ratio of the number of bridging siloxane (SiO-Si) bonds actually formed at a given time to the maximum number of such bonds expected for a completed reaction, that is, twice the number of silicon atoms present in the system. Parts a and b of Figure 1 show the variation of the degree of reaction with system density at low and high values of ω (1 and 5 respectively). For ω ) 1 (part a of Figure 1), we observe that over the entire range of system densities explored (0.2 to 1.8 g/cm3), the degree of reaction attains a fairly high equilibrium value of 80-90%. The slowest reaction kinetics was observed for the lowest-density configurations. These values of the degree of reaction are quite close to results from NMR studies of sol-gel polymerization of TEOS at 21 °C, where the degree of condensation was found to saturate at about 84%.13 Similar trends are observed for the ω ) 0 case. However, as the water-to-silicon ratio is gradually raised from ω ) 1 to ω ) 5, there is a considerable drop in the equilibrium degree of
1766 J. Phys. Chem. C, Vol. 112, No. 6, 2008
Figure 1. (a) Degree of reaction for bridging siloxane (Si-O-Si) bonds, ω (water-to-silicon ratio) ) 1 (densities in legend refer to overall system density, lines serve as guides to the eye); (b) degree of reaction for bridging siloxane (Si-O-Si) bonds; ω ) 5.
reaction. This drop is especially noticeable at relatively high system densities, where the degree of reaction stabilizes at about 50% or less (part b of Figure 1, ω ) 5). Thus, a high waterto-silicon ratio appears to deter the condensation reaction (1) at our simulation temperature (700 K). Most likely, the presence of additional water molecules in the system allows for differently structured hydration layers surrounding silicon-bearing species, which results in a shift of the equilibrium in eq 1 to the left encouraging the reverse (hydrolysis) reaction. Whereas this observation contradicts previous simulations carried out at 2500 K,41 29Si NMR studies of acid-catalyzed TEOS polymerization at room temperature actually found the gelation time to go up from 35 to 60 days as the TEOS/water molar ratio was reduced from 1:4 to 1:10.6 Experimental studies clearly show that a high water concentration promotes the hydrolysis of alkoxide precursors,6 but they typically do not reveal the effect of the water concentration on the rate of condensation of fully hydrolyzed precursor molecules, making a direct comparison with our simulation difficult. Some experimental studies1,55 indicate that the gelation time falls with increasing water-to-silicon ratio up to ω values of 4 to 6, beyond which it rises with increasing ω. This suggests that, at low ω values, the predominant effect of the addition of water may be to facilitate hydrolysis, whereas
Bhattacharya and Kieffer at higher ω values the addition of water mainly has the effect of slowing down condensation. III.1. Emergence of Distinct Growth Regimes. We observe a strong dependence of the structures that form on the system density and the water-to-silicon ratio. We can distinguish between three different growth regimes, characterized by unique structural features: (a) A small number of relatively large compact silica aggregates, or nanoparticles, prevail in systems with low densities for all water-to-silicon ratios. (b) Networked silica structures with percolation in one or more dimensions are observed at intermediate to high system densities and low water-to-silicon ratios. (c) A large number of branched silica clusters of varying size in solution with a low degree of reaction develop at intermediate to high system densities and high water-to-silicon ratios. A qualitative and quantitative assessment of the structures resulting from these growth regimes is given in parts a-c of Figure 2, which show representative examples of the three different structures along with the corresponding radial distribution functions (RDF) in the inset. Only silicon and oxygen atoms that have participated in condensation reactions are shown for clarity; line segments indicate Si-O bonds - the figures thus provide a visual impression of all of the silica aggregates present in the system. The configurations shown reveal the structures after the condensation reaction has reached equilibrium. For all of the growth regimes, the onset of the condensation reaction is characterized by the formation of a large number of small clusters, originating at various nucleation sites evenly distributed in space. In regimes (a) and (b), these clusters impinge upon each other as they grow, and fuse to form a small number of larger structures. Interestingly, the coalescence of particles in this fashion leaves little evidence of prior boundaries, that is, the silica networks merge uniformly, without encapsulation of condensed water molecules. In the example shown for regime (a), growth subsides with two distinct particles (part a of Figure 2), one containing 388 silicon atoms and the other, 121 silicon atoms. These two particles include almost all of the 512 silicon atoms present in the system. (Note that particles extend beyond the surfaces of the simulation box and reemerge on the opposite sides due to periodic boundary conditions; they therefore appear as two or more clusters in spite of being connected by chemical bonds.) This configuration is characterized by an overall density of 0.2 g/cm3 and ω ) 0. In regime (b), all of the initial growth clusters merge to form one continuous network. Part b of Figure 2 shows an equilibrium configuration resulting from this growth regime with a system density of 1.0 g/cm3 and ω ) 0. Apparently, the close spatial proximity of individual clusters during growth at higher density enables the continued aggregation. In contrast, at lower densities, once individual clusters have grown to a certain size and their mobility is reduced accordingly, the separation between clusters becomes increasingly more difficult to overcome, and aggregation comes to a halt on the time scales studied (i.e., ∼1 ns at 700 K). Part c of Figure 2 shows a typical example of regime (c): a large number of separate, branched clusters of varying size - in a configuration with a system density of 1.0 g/cm3 and ω ) 5. The differences in the structures resulting from each of these three growth regimes can also be recognized by comparing the respective Si-Si RDF (parts a-c of Figures 2, inset). Beyond the first- and second-nearest neighbor peaks, the Si-Si RDF corresponding to the configuration with compact silica particles shows a clearly defined maximum at about 1.5 nm, followed by a distinct minimum at about 2.3 nm (part a of Figure 2),
Polycondensation of Silicic Acid
J. Phys. Chem. C, Vol. 112, No. 6, 2008 1767 TABLE 1: Comparison of Average (Weighted) Radius of Gyration of Compact Silica Particles Observed in Systems of Low Density (0.2 g/cm3) with the Location of the Maximum in the Corresponding Si-Si Radial Distribution Function (rdf) water-to-silicon location of maximum in average (weighted) radius ratio, r SisSi rdf (Å) of gyration (Å) 0 1 2 3 4 5
Figure 2. (a) Configuration with compact silica particles (distinct particles colored differently, Si-Si radial distribution function shown in inset); (b) configuration with silica network structure; (c) configuration with branched silica clusters; (d) cumulative radial distribution function of structures corresponding to the three observed growth regimes, compared to that of simulated silica glass at 2.2 g/cm3.
15.2 11.4 9.1 13.2 9.1 11.5
15.7 11.9 9.1 10.8 9.4 10.9
indicating that particles are distinguishable and revealing their characteristic length scale. In case of the continuous network structure (part b of Figure 2) or the branched silica clusters (part c of Figure 2), the Si-Si RDF reveals no such characteristic length scale. For structures grown in regime (b), a weak undulation in the RDF hints at the fact that the continuous network arose by fusing together individual clusters, whereas in growth regime (c) the resulting structure is homogeneous as is evident from the uniformly concave shape of the RDF. The distinction between particles of regime (a) and clusters of regime (c) is also based on the difference in packing density between these aggregates. This is illustrated in part d of Figure 2, showing the cumulative Si-Si pair correlation functions, that is, the number of atoms found within a certain radius from any given atom, for structures corresponding to the three growth regimes as well as for dense silica glass as a benchmark. Accordingly, in the compact particles or the continuous networks, any given SiO4 tetrahedron is surrounded by other tetrahedra with a density comparable to that in silica glass up to at least the second nearest neighbor shell (∼0.5 nm), whereas in the branched clusters the packing density is much lower in this range. One measure of the characteristic size of the compact particles of regime (a) is the radius of gyration, RG, with respect to the center of mass of the particle, which can be directly calculated from simulation data. Whereas for structures grown in regime (b) the silica units assemble into one continuous network, structures grown in regimes (a) and (c) typically have between 2 and more than 10 separate aggregates. Radii of gyration of the particles from regime (a) ranged from 5.3 to 16.9 Å, depending on the reaction conditions. Some particles, grown under different conditions, had similar radii of gyration but contained very different numbers of atoms due to considerable variation in shape. To compare the information revealed by RDFs with the radii of gyration, we calculated an average of this latter measure, weighing the contribution of each particle by the number of silicon atoms it contained. We summarize the results in Table 1. Accordingly, the average radius of gyration closely matches the location of the local maximum in the corresponding Si-Si RDF. The only instance where this correlation is not very strong is in the system with ω ) 3; the possible reason being that there were only two particles present in this system and they varied greatly in size. We can thus conclude that the maximum in the Si-Si RDF is a good estimate of the average particle radius of gyration in systems made up of compact particles of comparable size. It is also worth noting that the average particle radius (Table 1) is similar to that observed in a recent study of silica nanoparticles formed in the presence of amino acids, specifically, by the addition of TEOS to an unbuffered aqueous solution of lysine.56 III.2. Structural Defects. As the condensation reaction progresses, the number of bridging Si-O-Si bonds in the system increases at the expense of Si-O-H bonds and
1768 J. Phys. Chem. C, Vol. 112, No. 6, 2008
Figure 3. Ratio (x) of the number of Si-O-Si bonds to the total number of oxygen atoms in the gel phase for different water-to-silicon ratios, ω, plotted against gel density (lines serve as guides to the eye).
unsaturated Si-O. The latter two bond types involve unreacted or partially reacted network oxygen. We denote by x the fraction of bridging oxygen relative to all of the network oxygen (i.e., not including the oxygen in water molecules). A high value of this ratio reflects a well-reacted system, whereas a low value indicates a high level of structural defects in the system. In particular, this quantity reveals what fraction of oxygen atoms constitutes the bulk versus the surface of a cluster or network fragment. Figure 3 shows x plotted against the effective density of the gel phase (i.e., the density calculated considering only the atoms making up the silica network) for different waterto-silicon ratios. The data points in Figure 3 appear to be scattered roughly into three different groups, which correspond to the three growth regimes (a), (b), and (c), described above. The first group, which comprises the low-density (
