J. Phys. Chem. B 2006, 110, 319-333
319
Generation of Atomistic Models of Periodic Mesoporous Silica by Kinetic Monte Carlo Simulation of the Synthesis of the Material Christian Schumacher,† Jorge Gonzalez,‡ Paul A. Wright,‡ and Nigel A. Seaton*,† Institute for Materials and Processes, School of Engineering and Electronics, UniVersity of Edinburgh, King’s Buildings, Mayfield Road, Edinburgh, EH9 3JL, United Kingdom, and School of Chemistry, UniVersity of St. Andrews, St. Andrews, Fife, KY16 9ST, United Kingdom ReceiVed: September 14, 2005; In Final Form: October 30, 2005
We have developed a molecular simulation method for the generation of realistic atomic-level models for periodic mesoporous silicas. Using simplified interaction potentials and simplified representations of the templating micelles, the simulation follows the reaction path of the hydrothermal synthesis and calcination of the silica material in a kinetic Monte Carlo (kMC) simulation. The only input to the simulation is the geometry of the micelle and the number of silicic acid monomers at the beginning of the synthesis. We simulated the adsorption properties of the PMS models using Grand Canonical Monte Carlo simulation. With use of MCM41 materials of different pore sizes as a prototype material, experimental and simulated adsorption isotherms for nitrogen, ethane, and carbon dioxide were compared, showing good agreement between simulation and experiment.
Introduction Recent progress in the synthesis of periodic mesoporous silica (PMS) materials enables the production of materials of various pore sizes and shapes in the mesoporous size range. PMSs are, at an atomic level, either amorphous or at most only partially crystalline, but nevertheless have a well-defined pore structure with long-range order.1-3 Our interest in PMSs is as adsorbents for gas separation and storage in industrial and environmental applications, and in particular in the design of these materials using molecular simulation-based tools. One element of the design process is the simulation of adsorption in model materials, usually using Grand Canonical Monte Carlo (GCMC) simulation.4,5 This might be done with one of two objectives: (i) predicting the performance of an existing adsorbent in new applications, or under changed conditions (thus minimizing the need for experimentation), or (ii) identifying an optimal structure for a particular application, which might involve choosing among different types of PMSs, different synthesis conditions, and the type of surface functional groups (if present). To achieve this, atomistic models of PMSs are required that are sufficiently realistic to give quantitatively accurate predictions of adsorption. Unlike in the case of, e.g., zeolites, the PMS structure is amorphous so a suitable model is one that is representative from a statistical point of view, rather than correct atom by atom. The realism of the model, in this context, is defined primarily by its ability to predict adsorption, though it should also be as consistent as possible with other information (such as XRD data and experimental data on the connectivity of the SiO2 network). The hypothesis we are testing in this work is that the most realistic models for PMSs are likely to come from a process that mimics the hydrothermal synthesis process used to make the real materials. We use MCM-41 as a prototype material in * Address correspondence to this author. E-mail:
[email protected]. Phone: (+44) 131 650 4867. † University of Edinburgh. ‡ University of St. Andrews.
this study, as it has a simple structure, and it has been widely studied. The surfactant-templated synthesis of periodic mesoporous silica is commonly carried out under acidic or basic conditions at ambient pressure and temperature. The mother liquor is an aqueous solution of a liquid-crystal forming surfactant, a source of silicic acid, and a pH-controlling acid or base. The periodic silica structure is formed by condensation of the silicic acid monomers to give a periodic silica skeleton that encloses the surfactant micelles. The morphology of the solid phase strongly depends on the composition, temperature, and pH of the precursor solution.6 The solid is obtained from the mother liquor by filtration and, to make the pore space accessible, the solid phase is heated to degrade the organic templates. Additionally, the high temperature during this last calcination step helps to increase the degree of polymerization (DP). Several variations of this basic synthesis procedure have been developed to obtain different geometries,7-11 improved long-range order, higher DP, and better hydrothermal stability.12-18 A generally accepted mechanism for the formation of the periodic mesophase during the hydrothermal synthesis of MCM41 is as follows.19 The surfactant molecules form rodlike micelles that are dispersed in the solution. Silicic acid monomers and oligomers interact with the micelles and through condensation a thin layer of silica is formed on the surface of the micelles.20 Subsequently, several of these rods aggregate in a regular array, which generates the long-range order observed for MCM-41. Further condensation cross links the aggregated rods and a stable periodic silica skeleton is obtained. The polymerization of silicic acid in different solutions has been studied by various experimental methods. 1H NMR and 29Si NMR experiments yield the distribution of silicon atoms of different numbers of siloxane bridges. The five different silica species are designated as Q0 to Q4, depending on the number of bridging oxygen atoms (bOs) and nonbridging oxygen atoms (nbOs) that are bound to the silicon atom. A Q0-silicon atom has no bOs but four nbOs while the opposite is true for a Q4-
10.1021/jp0551871 CCC: $33.50 © 2006 American Chemical Society Published on Web 12/09/2005
320 J. Phys. Chem. B, Vol. 110, No. 1, 2006 silicon atom. The distribution of the Qn-species gives evidence about the DP and the kinetics of the gelation reaction.21-24 Transmission electron microscopy, XRD, and small-angle neutron scattering experiments3 reveal the structure of the silica skeleton. Adsorption experiments with different adsorptives give information about the porosity, the pore size distribution,25 the surface area,11,26 and surface chemistry.27,28 The generation of atomistic models of PMS by a standard molecular dynamics (MD) simulation of the synthesis is not feasible because of the size of the system and the time scale that has to be covered by the simulation. Thousands of atoms are involved in the smallest representative section of the material, and the aggregation toward a regular phase takes minutes (and the calcination several hours), while MD simulations are in practice limited to the nanosecond time scale. Several approaches have been taken to circumvent this problem. For example, MCM-41 has been modeled as isolated one-dimensional channels with regular and disordered walls,29,30 as oxygen atoms distributed randomly around the pore volume,31 as pores cut out of solid amorphous silica obtained from the MD simulation of the annealing of molten glass,32 or by relaxation of randomly generated amorphous silica walls by MD simulation at high temperature.33 However, these models are not based on the physics of the reaction, so their ability to predict adsorption (without the tuning of the structure for particular adsorption systems) is likely to be limited. The polymerization of silicic acid monomers in the absence of templating micelles has been simulated in MD runs at high temperatures, where the reaction rate is sufficiently fast for the reaction to be accessible by MD,34-37 and in dynamic Monte Carlo simulations that do not include steric influences during the polymerization.38,39 Wu and Deem studied the nucleation of silica clusters from acid monomers in a Monte Carlo scheme that takes the atomic positions of the atoms into account.40 In addition, silica clusters and their interaction with other molecules present in the synthesis solution have been studied with ab initio methods.41-45 Larger systems of amorphous silica and silica surfaces have been studied over time spans of nanoseconds in MD simulations.28,46-49 The simulation time can be increased for small systems by combining Monte Carlo and Molecular dynamics techniques in a hybrid Monte Carlo simulation.50 Tu and Tersoff simulated amorphous silica as a continuous random network (CRN) in Monte Carlo simulations.51 The picture of a CRN is also valid for surfaces of amorphous silica.52 This approach has been modified to include nonbonded interaction between atoms and network defects such as undercoordinated silicon and oxygen atoms, to simulate the chemical vapor deposition of silicon and oxygen on a substrate.53 In this paper we present a kinetic Monte Carlo (kMC) method for the simulation of the hydrothermal synthesis and calcination of PMS. It enables the simulation, at an atomic level, of the entire process of the synthesis of templated PMS, starting from silicic acid monomers and the surfactant micelle and ending with the calcined structure, at the same temperature as the real synthesis process. As will be seen, our method depends for its practicability on some drastic (though physically based) approximations. As we will demonstrate, our approach allows us to make very accurate predictions of adsorption, and gives other material properties that are broadly in agreement with experiment, suggesting that our algorithm gives structures that are sufficiently realistic for this purpose. Experimental Section Synthesis of Calcined MCM-41. The reactants for the synthesis of MCM-41 were tetraethyl orthosilicate (TEOS, 98%
Schumacher et al. from Aldrich), hexadecyltrimethylammonium bromide (CTMABr, Aldrich) as a surfactant, sodium hydroxide (Fischer Chemicals), glacial acetic acid (Aldrich), and water. The molar compositions of the reactants used in these syntheses were 1:0.7: 0.12:589 for Si:NaOH:CTMABr:H2O, respectively. CTMABr was added to the aqueous solution of NaOH, followed by the addition of silica source, TEOS. After 20 h of stirring at room temperature the pH value was adjusted to 8.5 with acetic acid and the reaction medium was placed in an autoclave at 100 °C. After 24 h the solid was filtrated and dried at 60 °C in air for a further 24 h. The surfactant was removed by heating the solid at 550 °C under flowing nitrogen for 5 h and then in oxygen for another 5 h. N2 Adsorption Experiments. Nitrogen adsorption isotherms were measured at 77 K with an IGA-2 series Intelligent Gravimetric Analyzer. Samples were degassed at 110 °C overnight prior to measuring the isotherm. Adsorption values were measured gravimetrically. Each step was permitted to approach equilibrium over a period of 90 min. The samples were inserted untreated, i.e., the powder as obtained from the calcination. Ethane and CO2 Adsorption Experiments. The adsorption equilibrium isotherms for ethane and carbon dioxide on MCM41 were measured volumetrically at 263 K and up to 2.5 MPa with use of a bench-scale adsorption/desorption apparatus fitted with two Baratron absolute pressure transducers (MKS type 127A) with a two-channel readout/signal conditioner (MKS type PR4000). The reading accuracy is 0.05% of the usable measurement range (0-133.33 kPa and 0-3.3 MPa). The amount of sample used was about 4 g. The samples have been outgassed at 473 K for 3 h prior to the experiment. The ethane and carbon dioxide gases were supplied by BOC with a purity of 99.99%. Full details of the operating procedure can be found in ref 29. Computational Methods: Kinetic Monte Carlo The kMC method presented in this paper is an off-lattice simulation that follows the pathway of silica sol-gel processes in the presence of templating micelles. Covalent interactions between the silicon and oxygen atoms in the silica structure are treated by using the potential developed by Tu and Tersoff51 for a CRN description of bulk silica. This approach has been adapted to include nbOs and noncovalent atomic interactions, similar to the approach used by Burlakov et al. to simulate chemical vapor deposition of silica.53 Additionally, to be able to simulate the synthesis of PMS materials, we have allowed for the condensation reaction of nbOs and we have introduced models of the templating surfactant micelles. Our simulation focuses on the silica network, and ignores interactions that do not have a substantial impact on the formation of the network. Thus, the energy of the system is calculated from the positions and bonding of the bonded silicon and oxygen atoms, and the interaction with the micelles. This enables the generation of a realistic sequence of condensation reactions resulting in atomic models of a PMS structure, without simulating atomic interactions and motions that are not related directly to changes in the network topology. The advantage of this approach is that slow chemical reactions such as the polymerization of silicic acid at ambient temperature can be simulated over time scales far beyond those of MD simulations. Model for the Simulated System. The simulation volume is defined as a parallelepiped with periodic boundary conditions in all three directions as depicted in Figure 1a. The lengths of the edges are varied to minimize strain in the simulation. The
Simulation of the Synthesis of Periodic Mesoporous Silica
Figure 1. Schematic drawing of simulation cells with the defining vectors and angles (a) for an arbitrary parallelepiped and (b) for MCM41 including a representation of the micelle.
angles between the edges are arbitrary but they remain constant during the simulation. For a cubic simulation cell, all three angles are 90°. For the hexagonal lattice of MCM-41 one angle is 120° and two angles are 90°, as indicated in Figure 1b. The interaction between silica species, surfactant molecules, and other compounds in the real synthesis solution is very complex. In principle, an atomic representation of the surfactant chains is required to capture all relevant details. For example, water and silica are expected to penetrate the micelle more than 1 nm deep.54 However, to avoid the computational costs of such a detailed approach, a very basic description of the interaction between silica species and micelles is applied. The surfactant micelles, including possible counterions, which are templating the pores in the hydrothermal synthesis process, are defined in the simulations as simple geometric bodies. They interact with the oxygen atoms of the surrounding silica by an arbitrary, continuous potential function that is attractive for atoms that are outside the micelles, repulsive for atoms that penetrate the micelle, and zero for atoms that are further away than a specified cutoff radius. Assuming short diffusion pathways and small diffusional resistance to water in the relatively open network during both the hydrothermal synthesis and calcination, water molecules that are produced or consumed in the reactions are not represented explicitly. Furthermore hydrogen atoms are not taken into account directly, so the simulation does not distinguish between protonated and deprotonated forms of silanol groups. In the real synthesis, the pH value of the solution is adjusted by using hydrochloric acid or ammonium hydroxide, which is necessary to catalyze the hydrolysis of TEOS to silicic acid at ambient temperatures.34 The effect of the catalyst is sufficient to achieve a hydrolysis rate that is roughly 30 times faster than the rate of condensation of silicic acid monomers. Thus the fraction of unhydrolyzed TEOS in the synthesis solutions is negligible21 and can be ignored in the simulation. In addition to the catalytic function, the pH value influences the mechanism of the polymerization and acid-catalyzed systems show more open, less cross-linked colloidal particles than basecatalyzed systems that form more regular structures.26,34 This has to be taken into account implicitly in the reaction probabilities during the kMC simulation. Following the approach of a first shell substitution effect, the probability of a silicon being involved in a condensation reaction depends on its number of siloxane bridges.38,39 Interestingly, the kMC simulations show good agreement with experimentally observed kinetics for an acidic system21-23 if the reaction probability is the same for all nbOs, i.e., the reaction probability of a silicon atom is proportional to the number of its nbOs that can react with other nbOs nearby. This will be discussed in more detail in the results section of this paper. The kMC Algorithm. In the kMC simulation, four different types of reaction are carried out: the formation of a siloxane
J. Phys. Chem. B, Vol. 110, No. 1, 2006 321 bridge; the hydrolysis of a siloxane bridge; the swapping of a silanol group and a siloxane bridge; and the swapping of two siloxane bridges belonging to two neighboring silicon atoms. In addition, two other types of change take place: the random displacement of all atoms in the simulation cell and the insertion of an additional monomer into the simulation volume. Each step in the simulation involves one of these six “trials”. All trials are attempted with equal probability. Silicon atoms always remain fully coordinated with four oxygen atoms, but oxygen atoms are bound to one silicon atom if it is a nbO or two silicon atoms if it is a bO. First the type of trial is chosen randomly. Then the trial configuration is generated and its energy is computed by relaxing the positions of the atoms. Finally the change in energy compared to the original relaxed configuration of the network is calculated, and the new trial configuration is accepted or rejected depending on the difference in energy. The acceptance probability A(ofn) for a transition of the system from the old configuration o to the new configuration n is given by
[ (
A(ofn) ) min 1, exp -
)]
U(n) - U(o) kBT
(1)
The exponential term is the Boltzmann factor of the difference in energy between o and n at temperature T. If the detailed balance is obeyed, i.e., if the probability of attempting the reverse trial n f o equals the probability of attempting the original trial o f n, then this algorithm generates a Markov chain of random configurations with a Boltzmann distribution of the energy. In contrast to standard applications of the Metropolis algorithm in Monte Carlo simulation, we do not intend to calculate average properties of the equilibrated system but rather we use the Metropolis method to generate a random reaction pathway that avoids configurations of unrealistically high energy and leads to a metastable state confined by energy barriers (rather than quartz, which is the equilibrium structure of silica under ambient conditions). This implies that a kMC run samples only the part of all possible configurations that is near the reaction pathway, thus, unlike a conventional, equilibrium MC simulation, a kMC simulation is not ergodic.5 It is also the case that the kMC simulation does not fully satisfy the detailed balance,5 as the minimization of the energy of the network that is carried out in every trial leads to configurations that are not unambiguously related to the configuration generated before the minimization. Despite this limitation, we will see later that the kMC simulation method generates realistic reaction kinetics and model silica materials. The vibrational energy is neglected in the kMC scheme, i.e., only the fully relaxed silica network is considered and the temperature applied in the simulation is a “network temperature” TN.51 Since the distribution of the total energy over the kinetic energy of vibration and the potential energy of the network topology is unknown, TN can only be estimated, rather than calculated exactly. This has been done on the basis of test runs which revealed the effect of the temperature on kMC simulation. Below 300 K a relatively high DP is maintained and the structural changes in highly condensed silica are minor. Above ca. 500 K the energetic barriers to reorganization of the network topology are passable and at temperatures higher than 800 K the DP drops significantly. Consequently we chose 200 and 600 K as the network temperatures corresponding to ambient and calcination temperatures, respectively. To form a siloxane bridge, two reacting nbOs that are bound to different silicon atoms and are within a distance of 0.4 nm
322 J. Phys. Chem. B, Vol. 110, No. 1, 2006 of each other are chosen randomly. The formation of 2-fold rings is not allowed since the high strain44 makes them energetically unfavorable and they rarely exist in real silica, if at all.48,52,55,56 The cutoff of 0.4 nm was applied to yield a reasonable acceptance probability without affecting the detailed balance with the hydrolysis trial. That is, the separation of the nbOs generated in hydrolysis trials must be smaller than the cutoff radius applied in the condensation trial. This is the case in more than 99.9% of the hydrolysis trials. The two nbOs are replaced by one bO that binds to both silicon atoms, increasing their Qn values by one. The energy of the formation of a siloxane bridge from two neighboring silanol groups was set to -13.4 kJ/mol, according to quantum mechanical calculations for the reaction of two silicic acid monomers to a dimer in aqueous environment.57 This value is only an estimate and its effect on the outcome of the kMC simulations is influenced by several factors. The most important ones are the uncertainty regarding the network temperature and the difference in strains in small rings obtained with the potential for interatomic interaction used in this work (see below) and values obtained by ab initio calculations for small clusters.44 In simulations in which a micelle is present, the value of the energy-well depth of the micelle-oxygen interaction is added to the reaction energy to balance the fact that the produced water molecule (which is not modeled explicitly in the kMC simulation) would experience the interaction with the micelle instead of the eliminated nbO. The reverse trial to condensation is to attempt to hydrolyze a siloxane bridge. A bO is chosen randomly and removed. Two new nbOs reestablish the 4-fold coordination of the silicon atoms involved, reducing their Qn values by one. During the simulation of the hydrothermal synthesis, the reaction is reversible, thus condensation and hydrolysis are attempted with equal probabilities. In the subsequent simulation of the calcination at high temperature, the hydrolysis trial is disabled to reflect the lack of water present during the calcination process. The third trial, which changes the network topology, is the switching of two O-Si bonds of two silicon atoms that are connected by a third oxygen atom.58,59 A sequence O(1)-Si(1)O(2)-Si(2)-O(3) changes to O(3)-Si(1)-O(2)-Si(2)-O(1). Quantum mechanical simulations indicate that this is similar to the real diffusion mechanism in amorphous silica.60 To have an effect on the network structure, at least one of the switching oxygen atoms has to be a bO. The fact that the energy that goes into the calculation of the acceptance criterion is the minimized energy rather than the average of a short MD run limits the positions accessible to mobile small oligomers and flexible linear chains and they might remain in a metastable state for a long time. To reduce the effect of this shortcoming of the kMC method, a shake move was introduced, where all atoms are displaced in a random direction. The distance of the displacement was set to the length of the silicon-oxygen bond, 0.16 nm. Potentials. To be able to simulate the silica networks of mesoporous materials, the applied potentials have to be simple. The energy due to strain in covalent silicon-oxygen bonds is calculated by using the harmonic spring potential of Tu and Tersoff, which was fitted to reproduce the experimental bond length and bond angles and the elastic constant of amorphous silica.51 Noncovalent interactions are modeled by using a quadratic approximation of the repulsive part of the 12-6Lennard-Jones potential with parameters commonly used for oxygen atoms.
Schumacher et al. The templating micelles are represented by spheres or rods of the appropriate size, which attract nearby oxygen atoms, to mimic the electrostatic interaction between the silica and the charged surface of the micelles. Oxygen atoms that penetrate the template are repulsed, but the repulsion is soft enough to allow the surface of the micelle to adapt to small irregularities of the silica wall. The range of the attractive field is large enough so that all oxygen atoms are attracted by at least one template to counterbalance the absence of interatomic attractive forces. Silicon atoms do not interact with the templates since they are shielded by the surrounding covalent bound oxygen atoms. The total energy Utot consists of three parts as in eq 2. The energy due to stretching and bending of covalent bonds between silicon and oxygen atoms, Ubond, is calculated as in eq 3 by using the same parameters as Tu and Tersoff.51 The first sum runs over all Si-O bonds in the simulation cell, and the second over all bond angles between two bonds i and j of the same atom a. The equilibrium Si-O bond length b0 is 0.16 nm and the corresponding spring constant kstretch/kB is 31.33 × 106 K/nm2. The equilibrium bond angles θ0,O and θ0,Si are 180° and 109.47°, and the constants kbend,O/kB and kbend,Si/kB are 50131 and 8703 K, respectively. Note that the equilibrium Si-O-Si angle, θ0,O, differs from the maximum of the bond angle distribution obtained from experiments and quantum mechanical calculations which lie around 145°.60 If θ0,O had the experimental value, rings with three or more silicon atoms could relax any tension by rotations about the six Si-O bonds. With the larger value of 180 °, there is strain in any ring that has less than five silicon atoms, as is expected for real silica. The absolute values of the resulting energetic penalties are about 50% smaller than those derived from ab initio calculations for small cyclic silica clusters,44 but the relative differences obtained for three-, four-, and five-membered rings are the same. The unphysical Si-OSi bond angle will affect the orientation of the nbO atoms of small noncyclic oligomers, but since they are very mobile and since cyclization occurs early in the polymerization process, we assume that its effect on the final network topology is negligible.
Utot ) Ubond + Urep + Umic Ubond ) (1/2)
∑i kstretch(bi - b0)2 + (1/2)∑∑kbend,a(cosθij - cosθ0,a)2 a i,j
(2)
(3)
The repulsive term Urep for noncovalent interactions is calculated according to eq 4 for all pairs of atoms that are not bound to the same atom and within the cutoff radius rc of 0.28 nm. The constant krep/kB is 541.48 × 106 K/nm.4 The repulsive term does not affect the strain of rings of three or more silicon atoms. Since the formation of two-membered rings is not allowed in the kMC algorithm, the parameters for the covalent bonds control the energy penalty for strain in small rings.
Urep ) (1/2)
∑ krep(rc - rab)4
(4)
a*b rab
