Kinetic Monte Carlo Simulation of the Synthesis of Periodic

Aug 29, 2012 - *Phone: +44 (0) 1316504856; e-mail: [email protected]. ... was achieved demonstrating the realism of the pore models obtained by the ...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCC

Kinetic Monte Carlo Simulation of the Synthesis of Periodic Mesoporous Silicas SBA‑2 and STAC-1: Generation of Realistic Atomistic Models Carlos A. Ferreiro-Rangel,† Magdalena M. Lozinska,‡ Paul A. Wright,‡ Nigel A. Seaton,†,§ and Tina Düren*,† †

Institute for Materials and Processes, School of Engineering, The University of Edinburgh, The King’s Buildings, Mayfield Road, Edinburgh EH9 3JL, United Kingdom ‡ EaStCHEM School of Chemistry, University of St. Andrews, Purdie Building, North Haugh, St. Andrews KY16 9ST, United Kingdom S Supporting Information *

ABSTRACT: SBA-2 and STAC-1 are two related periodic mesoporous silicas (PMSs) that have regular networks of spherical, interconnected pores; the pores are similar in the two materials but the networks differ in their symmetry. The nature of the interconnected network of pores in these materials gives rise to interesting properties related to their potential use in separation processes. In this work, we extend a kinetic Monte Carlo (kMC) technique, originally derived for MCM-41, a simpler PMS, and apply it to mimic the condensation, aggregation, deformation, and calcination stages of the synthesis of SBA-2 and STAC-1. Our simulated synthesis results suggest that the pores are connected through windows formed during micelle aggregation because of the close packing of the spherical micelles and the presence of water molecules at the silica−micelle interface. The simulated materials were validated by comparing properties such as unit cell size, pore size, pore shape, and wall density to results from experimental X-ray diffraction (XRD), transmission electron microscopy (TEM), density measurements, and 29Si NMR. Quantitative agreement between simulated and experimental nitrogen isotherms was achieved demonstrating the realism of the pore models obtained by the kMC simulations. Our results highlight the importance of a realistic, rough pore surface for the prediction of adsorption at low pressures in these materials.

■

INTRODUCTION Periodic mesoporous silicas (PMSs) are porous materials that show long-range order but are amorphous on the molecular scale. They have been widely studied since 1990 when the Mobil Corporation Laboratories published the controlled synthesis of MCM-41.1 MCM-41 is a material with cylindrical pores arranged in a hexagonal array and remains one of the most widely studied PMSs.1−10 Since the first synthesis of MCM-41, many PMSs with different pore shapes and structures have been synthesized. There is a special interest in these materials since properties such as pore size, shape, and connectivity can be tailored by adjusting the synthesis conditions. Additionally, the pore surface of PMSs can be modified postsynthesis to optimize the material for specific applications.11−13 In this paper, we focus on SBA-2 and a material with a related structure known as STAC-1.14 SBA-2 was first synthesized by Huo et al. in 1995.15 It is a mesoporous material with spherical pores and P63/mmc symmetry resulting from the hexagonal close packing (hcp) of spherical surfactant micelles. The ratio of its unit cell parameters (c/a) is typically between 1.61 and 1.63.16 This © 2012 American Chemical Society

material can be represented by an infinite two-layer sequence A-B-A-B... of spherical pores in which each pore has a coordination number (i.e., the number of neighboring pores) of 12. SBA-2 is produced by using gemini quaternary ammonium surfactants. As for other PMSs, the length of the hydrophobic tail of the surfactant used has a direct effect on the pore size of the final material, which exhibits a very narrow pore size distribution.15 The nitrogen adsorption isotherm in this material is of IUPAC type IV showing that SBA-2 pores are in a connected network. Some work has been done to understand the nature of the pore network in these materials.17 It has been suggested that pore connectivity results from the organic compounds breaking their way out of the material during the calcination stage,18 but this remains an unproven hypothesis. Pérez-Mendoza et al.17 represented the SBA-2 structure as a system of spherical cavities with smooth (i.e., atomistic but regular) walls, and similarly Received: August 1, 2012 Revised: August 29, 2012 Published: August 29, 2012 20966

dx.doi.org/10.1021/jp307610a | J. Phys. Chem. C 2012, 116, 20966−20974

The Journal of Physical Chemistry C

Article

smooth cylindrical connecting channels, for which they fitted a pore size distribution (PSD) using experimental ethane isotherms. The PSD found was bimodal with two peaks at ∼9 Å and ∼47 Å representative of the channels and the cavities, respectively. An important limitation of their approach is the smoothness of the pore models used to represent SBA-2 since it is well-known that the pore surface of real PMS materials is rough.10,19 A complexity arises when synthesizing SBA-2, which is that it is also possible to obtain coexisting regions of hexagonal and cubic (Fm3m) symmetry as presented by Zhou et al.14 This latter symmetry leads to a slightly different structure named STAC-1, which is defined by a stacking sequence of three layers A-B-C-A-B-C... consistent with that of a face-centered cubic structure, that is, cubic close packing (ccp) of the spherical cavities. In both cases, the pore-packing efficiency is the same as is the coordination number (12). It is the presence of the third layer (C) that changes the structure by introducing an additional plane of symmetry obtained by rotating the structure about the vertical space diagonal of a cubic unit cell. Our aim is to carry out an atomistic simulation of the formation of SBA-2 and STAC-1 both as an aid to understanding the synthesis process and to create realistic models for the prediction of adsorption in these materials. In previous work, the kinetic Monte Carlo (kMC) technique was used to generate atomistic models for MCM-41 that were successfully used to predict quantitatively the adsorption of various gases.20−22 Here, we extend this kMC technique to create realistic models of the more complex SBA-2 and STAC-1 materials. This allows the study of the different pore connectivity in the SBA-2 and STAC-1 regions, which was first observed by Zhou et al.14 The model materials are validated by comparing their properties (such as unit cell size, pore size, pore shape, and wall density) to results from experimental X-ray diffraction (XRD), transmission electron microscopy (TEM), density measurements, and 29Si NMR as well as by predicting nitrogen adsorption at 77 K and comparing this isotherm to experimental measurements.

surfactant−monomer interactions. While for MCM-41 the micelle is a cylindrical rod,20 for SBA-2 and STAC-1, a spherical micelle is required. The focus of the kMC technique is the polymerization reaction that leads to the formation of the material rather than the physical interactions taking place between all atoms. The nonbonded energy contributions in the system come mainly from oxygen and micelle interactions, while silicon interactions are not taken into account (because they are always surrounded by oxygen atoms, their contribution is assumed negligible). Furthermore, despite the polar nature of the atoms involved in the synthesis (silicon and oxygen), their electrostatic interactions are not taken into account. This may have an effect on the synthesis path followed to generate the model materials, but it is unlikely to affect the realism of the final configuration as shown by the MCM-41 model materials obtained by Schumacher et al.20 The steps used in the kMC simulation reflect the microscopic events that occur during the real synthesis. These are condensation (formation of silica bonds), hydrolysis (breaking of silica bonds), swapping of bonds (between a bridging oxygen and a nonbridging oxygen of two neighboring silica monomers), and shaking (random displacement of the atoms in the unit cell). All these trials take place at all stages in the synthesis except for hydrolysis, which does not occur during calcination.20 Simulations are set up by imposing a unit-cell shape and size as well as by specifying the number, size, and location of the micelles and the number of silica monomers in the simulation cell. The silica monomers are then placed outside the volumes occupied by the micelles. The first stage comprises the formation of a silica monolayer around the micelles at a temperature of 300 K. (This temperature is actually a network temperature, i.e., the temperature of the fully relaxed silica network, since the vibrational energy is not taken into account by the kMC technique,20 and thus, it is not possible to know its kinetic energy contribution.) During this stage of the experimental synthesis, the micelles have not yet formed an ordered array, and so the simulated micelles are treated as independent of each other, and periodic boundary conditions are not used. Unlike the previous work on the simulated synthesis of MCM-41,20 the kMC simulation of the synthesis of SBA-2 and STAC-1 explicitly accounts for water molecules (in the form of unbound oxygens), since it was found that water molecules have a direct effect on the silica polymerization by promoting the reaction and thus play an important role especially in the early stages of the silica structure formation. Immediately after the formation of a monolayer, the micelles aggregate to form a regular array and in the process undergo deformation. In the simulation, periodic boundary conditions are used to allow micelles to interact with the neighboring micelles and the silica layers that surround them. The reaction that follows creates the silica skeleton for the material. By using the appropriate unit-cell shape, we are able to obtain a realistic representation of the symmetry of the different STAC-1 and SBA-2 structures (as will be illustrated in the next section). After the completion of this stage (i.e., the equilibration of the rate of successful condensation and hydrolysis trials), the calcination stages take place. In the real system, there is a progressive disintegration of the micelles upon calcination. As the kMC simulation does not include an explicit account of the disintegration of the micelles, this is approximated by using two

■

COMPUTATIONAL METHODOLOGY Experimental Synthesis of SBA-2 and Its Representation by kMC Simulation. We adopt the kMC method developed by Schumacher et al.20 Their original paper gives a full description of the kMC algorithm. In this paper, we outline the basic elements of the kMC method and give an account of the adaptation needed to simulate the synthesis of SBA-2 and STAC-1. The kMC technique is an off-lattice simulation that evolves stochastically following a set of rules that mimic the sol−gel synthesis of PMSs via the condensation and aggregation mechanism.23 The experimental synthesis of PMSs goes through several stages−monolayer formation, aggregation and deformation, calcination, and cooling−all of which are represented in the kMC simulation. The kMC simulations are carried out using a simplified interaction scheme, where the silicon atoms are always fully coordinated to four oxygen atoms (thus, the exact reaction path is not mimicked) and where the hydrogen atoms are not explicitly modeled. Furthermore, since the shape of the templating micelle is responsible for the final shape of the cavities in the material, a geometric shape representative of the precursor micelles observed in the experimental synthesis is used thus avoiding the time-consuming calculation of 20967

dx.doi.org/10.1021/jp307610a | J. Phys. Chem. C 2012, 116, 20966−20974

The Journal of Physical Chemistry C

Article

in the simulation. A rhombohedral unit cell (i.e., α = β = γ = 60°) with a single micelle in it leads to the ccp symmetry exhibited by the STAC-1 structure (see Supporting Information). While preserving its shape, the unit cell and its micelle deform because of the strain introduced by the formation of the silica network. Compared to the simulation setup for STAC-1, the setup for SBA-2 differs in three ways to achieve the hcp packing: the shape of the unit cell is hexagonal (α = β = 90° and γ = 60°), two model micelles are placed in it (see the Supporting Information), and the c/a cell ratio is kept as closely as possible to the ratio of the hexagonal closed packed systems (∼1.62). In other words, the generation of the simulated SBA-2 material requires the explicit representation of the A and B layers in the unit cell, while for the creation of the simulated STAC-1 material, it is sufficient to represent one layer with the others (B and C) being obtained as a result of the periodic representation of the unit cell. Characterization of the Simulated Materials. Solid state 29 Si MAS NMR (magic angle spinning nuclear magnetic resonance) provides information on the connectivity of the silicon atoms in a material. The same information can be determined for the simulated materials by scanning the connections of the silicon. This is done throughout the simulated synthesis. Silicon monomers are then classified according to the number of other monomers to which they are linked as Q0, Q1, Q2, Q3, or Q4 as polymerization evolves from isolated monomers (Q0) up to full coordination (Q4). The silica structure is composed of rings, which are defined as the shortest nonreversing path of oxygen−silicon bonds connecting two different bonds on the same silicon atom.20 The number of silicon atoms in the ring determines its size; this is also obtained throughout the simulated synthesis by scanning the silicon atoms. To calculate the size of the simulated pores, we used a method based on a random walk. A probe sphere is initially placed at the center of the pore from where it moves with random orientation in steps of a size short enough to ensure that the new position overlaps the old one (i.e., not larger than the probe radius). This process is repeated for a large number of steps (two million, in our case). When the probe touches a wall atom, its distance to the center of the pore is reported before it is replaced in its previous position. The touched wall atom is tagged to avoid repeating measurements. Thus, only wall atoms on the pore surface are used in the pore radius calculation. Because of the rough nature of the pore surfaces, this method leads to a pore diameter distribution. We calculated the fourth central moment (the kurtosis) of the pore diameter distribution as this gives an indication about the surface roughness. Any material with a kurtosis larger than three is considered to be very rough.25 The pore diameters for the simulated materials presented in this work were obtained with a spherical probe molecule of d = 3.30 Å consistent with the Lennard-Jones diameter of a nitrogen atom since this is a fluid commonly used in the experimental determination of porosity. We also calculate the accessible surface area of the pores rolling a nitrogen-sized probe molecule across the surface.26 Finally, we report the density of the siliceous skeleton of the simulated materials. This calculation requires the mass of the solid as well as the volume occupied by the silica network (which is the total volume of the solid minus that of the pore space). Mimicking the experimental procedure where the pore volume is determined from the adsorption of a nonadsorbing

calcination stages in succession: one with the micelles still present followed by another where they are absent. In the first calcination stage (in which the micelles are still present), the network temperature is increased by 500 K (i.e., from 300 to 800 K), and the water molecules are removed. This means that the hydrolysis trials can no longer take place as there is no water to allow the reverse condensation reaction. Once the percentage of Q3 and Q4 silica monomers stabilizes, the micelles are removed and the second calcination stage begins. The second calcination stage ends once the acceptance of the condensation and swapping trials is almost negligible (i.e., the silica structure no longer varies). A final cooling stage (i.e., lowering the network temperature to 300 K and allowing hydrolysis) is then simulated; the silica network relaxes mimicking the cooling of the real material to ambient conditions. Throughout the application of the kMC technique, the potential energy of the system is calculated for the acceptance/ rejection of the trials. To this end, a network energy minimization takes place after each trial through a steepest descent algorithm thus kinetically rearranging the atom positions. As in the work by Schumacher et al.,20 the contributions to the potential energy of the system come from the strain of the Si−O bonds (modeled through the harmonic spring potential of Tu and Tersoff24), the repulsion experienced by close oxygen atoms (modeled by a quadratic approximation of the repulsive part of the 12−6 Lennard-Jones potential), the interactions between the oxygen atoms and the geometric model micelles, and the repulsion between neighboring model micelles. The real micelles have flexible rather than geometrically fixed boundaries. Schumacher et al.20 represented this by using a soft potential for oxygen−micelle interactions, which we also used in this work. This allows monomers to penetrate the geometric shape of the model micelles, which as we show below has an important effect on the surface roughness. Unlike the simulated synthesis of MCM-41,20 the repulsive interaction between micelles in the SBA-2 and STAC-1 systems needed to be stronger (Figure 1) to avoid the merging of different micelles, which was particularly evident in the simulated synthesis of SBA-2. Generation of the Model Materials. Whether STAC-1 or SBA-2 is generated depends on the shape of the unit cell used

Figure 1. Comparison between the micelle−micelle repulsion energies for spherical micelles with a radius of 30 Å: this work (solid line), for cylindrical micelles for the simulation of the synthesis of MCM-4120 (dashed line). 20968

dx.doi.org/10.1021/jp307610a | J. Phys. Chem. C 2012, 116, 20966−20974

The Journal of Physical Chemistry C

Article

gas such as helium, we carried out grand canonical Monte Carlo (GCMC) simulations of helium adsorption in the simulated materials using the fluid critical properties27 and the interaction parameters given in Table 1. The pore volume (VHe) was calculated using the ideal gas equation.28 Table 1. Interaction Parameters Used in the GCMC Adsorption Simulations pore wall

He N2

site

εi/kb [K]

σi [Å]

qi [eo]

ref

Si bO nbO H He Q N N

0.0 185.0 185.0 0.0 10.9 0.0 34.897 34.897

0.0 2.708 3.000 0.0 2.640 0.0 3.3211 3.3211

1.2805 −0.6402 −0.5261 +0.2060

20, 21, 5

+1.0950 −0.5475 −0.5475

Figure 2. Degree of polymerization, Qn, as obtained for a STAC-1 model material with pore radius of ∼20.6 Å at the end of the individual synthesis stages.

31 32, 33

the silica structure. The increased proximity of the silica monomers leads to additional silica condensation thus promoting the formation of the silica skeleton of the structure. The strain caused by the latter influences the deformation of the unit cell and the micelle radius during this stage (as occurred in the simulated synthesis of MCM-4120). During deformation, there is typically an initial increase in the micelle radius that maximizes the silica−micelle interface (i.e., minimization of the system energy), which is eventually replaced by shrinking of both the unit cell and the micelle radius as the strain of the silica structure increases. By the end of this stage, polymerization is well under way as can be seen from the high percentage of Q3 and the absence of isolated monomers in Figure 2. The first calcination stage has a predominant effect on the silica polymerization as can be seen from Figure 2. By the end of the first calcination stage, nearly all silica monomers are almost equally likely to be in either the Q3 or Q4 states with the nearly total absence of both isolated and small-clustered silica monomers. Figure 2 shows that, in the STAC-1 model material, by the end of the second calcination stage around 60% of the silica monomers are fully polymerized (Q4) with the remaining monomers belonging to the Q3 population. Figure 2 also shows that the further relaxation of the system during the cooling phase makes little contribution to the final characteristics of the simulated material. This confirms the stability of the simulated material since during this stage both condensation and hydrolysis can take place reflecting the fact that cooling a sample takes place under normal atmospheric humidity. In Figure 3, we present the Qn distribution for the simulated as-prepared STAC-1 material (i.e., after micelle aggregation) and the same material after calcination and cooling along with the experimental degree of polymerization. The experimental Qn distributions are obtained by deconvolution of experimental 29 Si MAS NMR spectra (shown in Figure S2 of the SI) of SBA2 samples which are always mixed-phase materials consisting of both STAC-1 and SBA-2. There is a very good agreement between the results from the simulations and those obtained experimentally although it can be seen that the simulated material shows slightly lower Q4 and higher Q3 in both stages. In particular, the percentage of Q2 for the simulated asprepared material is higher than that reported experimentally. Thus, the silica monomers in the experimental sample are less branched than those in the simulation cell at this stage. This

Grand Canonical Monte Carlo Simulation of Adsorption. We use grand canonical Monte Carlo (GCMC) simulation to calculate the equilibrium adsorption properties of the model material. In this method, the volume, temperature, and chemical potential are fixed, and the number of fluid molecules varies stochastically. The mean number of molecules inside the pore as a function of pressure gives the simulated adsorption isotherm. The Peng−Robinson equation of state29,30 is used to relate the chemical potential to the pressure in the system with the critical properties of the fluid taken from the NIST database.27 The GCMC simulation then allows the number of fluid molecules in the system to fluctuate thus giving the absolute amount adsorbed in the system at the given thermodynamic conditions. The absolute amount adsorbed, Nabs, relates to the experimentally measured excess amount adsorbed, Nex, by31 Nex = Nabs − ρbulk VHe

(1)

where ρbulk is the density of the bulk gas calculated using the Peng−Robinson equation of state, and VHe is the pore volume determined from simulated helium adsorption isotherms. The dispersive interactions in the system were modeled by the Lennard-Jones (LJ) potential, where σ and ε are parameters corresponding to the LJ diameter and potential well depth, respectively.28 The mixed interactions between atoms of different species were obtained using the Lorentz−Berthelot combination rules,28 and the long-ranged Coulombic interactions were modeled using the Ewald summation technique.28 The interaction parameters for the mesoporous silicas are taken from previous works on MCM-41.5,20 The Lennard-Jones parameters and partial charges used in the GCMC simulations are summarized in Table 1.

■

RESULTS Figure 2 shows the degrees of polymerization at the end of each simulated synthesis stage for the simulated STAC-1 material. The results for the simulated SBA-2 material are similar and are shown in Figure S1 of the Supporting Information (SI). Figure 2 shows that the transition between the monolayer formation and the aggregation stages is characterized by a low percentage of isolated monomers (