{'fraction_non_alphanumeric': 0.02442822047680011, 'fraction_numerical': 0.007075879727345075, 'mean_word_length': 4.879666361509511, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Program Report COMP 1 Modeling spontaneous formation of precursor nanoparticles in clear-solution zeolite synthesis

We present a lattice model describing the formation of silica nanoparticles in the early stages of the clear-solution templated synthesis of silicalite-1 zeolite. Silica condensation/hydrolysis is modeled by a nearest-neighbor attraction, while the electrostatics are represented by an orientation-dependent, short-range interaction. Using this simplified model, we show excellent qualitative agreement with published experimental observations. The nanoparticles are identified as a metastable state, stabilized by electrostatic interactions between the negatively charged silica surface and a layer of organic cations. Nanoparticle size is controlled mainly by the solution pH, through nanoparticle surface charge. The size and concentration of the charge-balancing cation are found to have a negligible effect on nanoparticle size. Increasing the temperature allows for further particle growth by Ostwald ripening. We suggest that this mechanism may play a role in the growth of zeolite crystals.


I. INTRODUCTION
Zeolites are nanoporous alumino-silicates used in industrial applications such as catalysis and separations. 1 Understanding how zeolites nucleate and grow is of fundamental scientific and technological importance. Such understanding could be used to optimize catalysis and separations by tailoring zeolite crystallite size and shape. In addition, controlled zeolite synthesis could open new fields of application such as optical electronics, 2 bio-implants 3 and enantioselective separations. 4 In the past decade, silica nanoparticles have been found to play an important role in zeolite formation from clear solutions. [5][6][7][8][9][10] To explore the properties of such nanoparticles, we have developed and applied a lattice model to simulate nanoparticle formation, structure and stability.
The clear-solution synthesis of silicalite-1 zeolite provides an important test case.
Silicalite-1 is the pure-silica form of zeolite ZSM-5 (framework code MFI), which was first reported by Flanigen et al. 11 It is composed of parallel straight cylindrical channels, connected by cylindrical zigzag channels. 12 This unique structure enables numerous applications of MFI-type materials, particularly in hydrocarbon conversions. 13 Silicalite-1 can be synthesized from a clear solution of tetraethyl orthosilicate (TEOS), tetrapropylammonium hydroxide (TPAOH) and water. When this solution is heated to ∼100 o C for a few hours, it yields high-quality crystals of silicalite-1. 14 Directly after synthesis, its channel intersections are occupied by tetrapropylammonium cations (TPA), with the four propyl groups extending into the adjacent channels. Removal of TPA by calcination produces a nanoporous crystalline solid.
During the clear-solution synthesis of silicalite-1, a suspension of nanometer-sized silica particles has been observed by several researchers prior to the appearance of zeolite crystals. [5][6][7][8][9][10]15 These nanoparticles are found to remain stable at room temperature, and lead to silicalite-1 crystals at elevated temperatures. Experimentally measured growth rates and activation energies were found to agree with the predictions of a mathematical model assuming that the rate limiting step for silicalite-1 growth is the addition of nanoparticles. 15 Even though the presence of these particles is agreed upon, their structure and their precise role in zeolite synthesis remain the subject of debate. Schoeman suggested that the nucleation and growth of the zeolite proceeds by dissolution of the nanoparticles and subsequent monomer and oligomer addition, 5 while de Moor et al. proposed a mechanism of direct ag-gregation of particles to form zeolite crystals. 7 Martens and co-workers have recently argued that the particles are pre-formed nanometer-sized slabs of silicalite-1, and that silicalite-1 crystals are formed by organized assembly of these units. 8 However, subsequent experiments have thrown doubt on these conclusions, suggesting that the nanoparticles are amorphous entities. 9 More recently, Fedeyko et al. published a detailed in situ small-angle x-ray scattering (SAXS) and small-angle neutron scattering (SANS) study of silica nanoparticles formed in basic solutions of tetraalkylammonium cations (TAA). 16 They interpreted their results in terms of a core-shell structure for the particles, with a core of mostly amorphous silica and a shell of TAA mainly located on the particle surface. They found that average particle sizes decrease with pH, increase with temperature, and are nearly independent of the size of the alkyl group in TAA. A decrease in particle size with pH has also been reported by Yang et al. 10 Despite the relative simplicity of this zeolite growth system, e.g., no alumina source, we lack a fundamental atomistic understanding of the basic interactions that control the structure and size of these nanoparticles. The objective of this paper is to develop a statistical mechanical model of the spontaneous formation of silica nanoparticles from an aqueous solution of TEOS and TPAOH. For the most part, we focus on the solutions studied experimentally by Fedeyko et al., 16 because they present a systematic analysis of several physical parameters. expedient calculations.
There have been several previous attempts to model silica polymerization using molecular simulation methods. [17][18][19][20][21][22] However, most of these simulations were performed at extremely high temperature, to allow for chemical bond breaking and reforming. Although valuable, this approach is inadequate for the system studied here, which deals with polymerization of silica at room temperature. Other approaches have been used to circumvent this problem, mainly based on employing "smart" Monte Carlo (MC) moves to sample the three-dimensional connectivity of silica polymers. 21,22 However, the current state-of-the-art in terms of computer power restricts these simulations to relatively small systems and/or short run times. In order to model the spontaneous formation of nanoparticles of realistic sizes, and to obtain statistically meaningful information about their properties, one must sample over large system volumes during long simulation times. Such simulations are currently too computationally intensive for off-lattice atomistic models. Lattice models, on the other hand, are capable of providing at least qualitative information when applied to such complex systems. In particular, lattice models have been successful in predicting the phase behavior and structure of micellar solutions, [23][24][25][26][27] and have recently been extended to surfactant-silica systems. 28 Inspired by these successes, we have applied a lattice model to simulate nanoparticle formation during the clear-solution synthesis of silicalite-1.
The sol-gel chemistry of aqueous silica solutions exhibits a rich variety of behavior resulting from an interplay of several phenomena such as solvolysis, acid-base equilibrium, metastability and phase separation. [29][30][31] Furthermore, the high alkalinity of these solutions and the presence of organic cations means that electrostatics likely plays an important role.
The organic cations are also associated with a structure-directing effect in the synthesis of zeolites, and thus steric effects might also be important in the context of nanoparticle formation. In this paper, we boil down all this complexity into a simple lattice model. We achieve this by representing silica polymerization by an effective attraction between lattice sites, including orientation-dependent interactions meant to embody electrostatics, and using second-neighbor interactions to model size-exclusion effects. In addition to MC simulations in the canonical ensemble, we employ parallel tempering MC 32,33 to probe the existence of metastable states, as well as reactive ensemble MC 34,35 to study reaction equilibrium. We observe in our simulations the spontaneous formation of silica nanoparticles under conditions where they are seen experimentally. In our model, these particles are found to be a metastable state of the system, stabilized by the interactions between ionized silicates and organic cations. The model correctly reproduces several qualitative features of the experiments outlined above.
The paper is organized as follows. Section II describes the model we utilize for the TEOS-TPAOH-water system, and relates the model parameters to the experimental thermochemistry of silica solutions. Section III describes the simulation techniques employed and Section IV presents the results obtained from the simulations. The results section begins with a brief analysis of the parameter space of the model, and then describes the formation of silica nanoparticles over the course of a simulation. The conditions for the existence of this nanoparticle phase are established and the effect of several parameters on the particle size and distribution are determined. In Section V, we present our conclusions and a discussion of future work.

II. MODEL DESCRIPTION
Our goal is to describe nanoparticle formation using a simple statistical mechanical model. This is particularly important for nanoparticle formation, since large system sizes are necessary to obtain statistically meaningful information, and long simulations are required to observe the spontaneous formation and evolution of the particles. Another advantage of using a simple model is that the physical properties responsible for a given behavior of the model can be more easily identified. Obviously, we must pay a price for using such a coarsegrained model, and that price is the difficulty in obtaining quantitative insights about the system. The simplifications involved are such that one is usually restricted to qualitative comparisons to real systems. Keeping this in mind, we will endeavor to keep our model as simple as possible, while retaining a realistic enough representation of the physico-chemical processes at play. We thus begin by abandoning an atomistic representation in favor of a "united atom" representation. This means that a given molecule will be represented by a single site, rather than by an explicit representation of all its atoms. Furthermore, we divide three-dimensional space into a simple cubic lattice, with every lattice site occupied by a single species.
The starting point for the silicalite-1 synthesis is an aqueous solution of TEOS and TPAOH at room temperature. The control variables in the experimental system are therefore the temperature, the initial concentration of TEOS and the initial concentration of TPAOH (which sets the initial pH of the solution). The base dissociates completely into TPA cations and hydroxyl anions. 36 At the conditions commonly used for silicalite-1 synthesis (high water/silicon ratio and high pH), hydrolysis of TEOS, yielding silicic acid and ethanol, is complete. 30,31 We are thus left with a system containing water, ethanol, TPA, Si(OH) 4 and OH -. To simplify this picture, we begin by ignoring the effect of ethanol and consider the water/ethanol mixture as a single solvent species (hereafter designated by W).
In reality, the presence of ethanol may affect the silica polymerization equilibrium (by virtue of esterification reactions) and alter the solubility of silica (by changing the nature of the solvent). 29 Nevertheless, we assume that in the dilute solutions studied here, these effects are of minor importance. In addition, we do not explicitly represent hydroxyl ions, but rather include their effects in acid-base reactions described below.

A. Neutral Polymerization
In order to build our model, we begin by considering the simplest case possible and proceed by adding subsequent levels of complexity, until we reach a satisfactory representation of the real system. The simplest relevant case one can envisage is pure silica immersed in solvent at the isoelectric pH of silica. In this situation, the relevant chemical reaction (again, ignoring the presence of ethanol) is silica condensation/hydrolysis (here written in a generic form): A solution that is supersaturated with silicic acid monomer will quickly evolve into a complex solution containing multiple silicate species, including dimers, trimers, rings and branched structures. 29 We treat this mixture by beginning with neutral silica monomers (denoted by S N ), each occupying one site on the lattice. The polymerization reaction energy is mimicked by introducing an effective nearest-neighbor interaction between silica monomers.
For example, if during the simulation two monomers move close to each other and become separated by only one lattice bond, we say that a dimer has been formed. Similarly, if a third monomer moves to within one bond of any of the previous two, a trimer is formed, and so on. This move is attributed a certain energy, which represents the internal energy of the condensation reaction (1). This should be distinguished from the much larger energy of forming a covalent bond, i.e., a covalent association energy. Both computational 37 and experimental 29 studies have determined the silica condensation reaction to be exothermic, which means that the interaction energy in our model should be attractive.
Using this procedure, a system composed of only neutral silica species immersed in solvent is fully specified by a single energy parameter -the value of the attraction between silica monomers, relative to the system temperature. To see that this is so, we consider the Hamiltonian of such a system in the canonical ensemble by fixing the temperature T , the number of lattice sites L 3 (essentially volume), and the numbers of S N and W molecules.
For simplicity, we include only nearest-neighbor interactions, but the approach can be generalized for longer-range interactions. The Hamiltonian of such a system can be expressed as: where N ij is the total number of nearest-neighbor contacts between components i and j, ε ij is the energy of interaction between those species, and c is the total number of components in the system. It can be shown 28 that, when all sites are occupied by either solute or solvent, the energy difference between two configurations is given by: This means that the system is fully specified by a set of interchange energies ω ij , with i = j.
For the simple system composed of silica and solvent molecules, we need only specify a single interchange energy. If we now consider the case where all interactions involving the solvent molecules are set to zero, hence defining the zero of energy, we can see that the interchange energy is given simply by: Thus, the behavior of the system depends only on the value of the silica-silica interaction energy. Such a parameter should be interpreted as an effective interaction, lumping together not only the reaction energy for silica polymerization, but also the solvation energies of the participating species.
We calibrate this effective interaction by simulating the solubility of amorphous silica at its isoelectric point, as described in section IV A. In the remainder of this paper, all energies are scaled by |ε S N S N | (the absolute value is used because ε S N S N is negative and we wish to retain the sign of the interaction). The reduced temperature is given by where k B is Boltzmann's constant.
A corollary of equation (3) is that a given solute-solvent system, with a given set of interactions, can always be reduced to an equivalent system where all the interactions involving solvent are reduced to zero. In practice, this allows us to treat the sites occupied by solvent molecules as vacant sites. 38 As we will see in section III, treating the solvent sites as vacancies also allows one to speed up the calculations significantly.

B. Alkaline Polymerization
The procedure described above yields a simplified description of the equilibrium represented by reaction (1) under isoelectric conditions. However, at high pH an additional reaction becomes important: To include this reaction, one must distinguish between neutral and ionized monomers. A further simplification introduced in our model is to consider only singly ionized monomers, i.e., Si(OH) 3 O − . Even though doubly ionized monomers (Si(OH) 2 O −2 2 ) have been shown to exist at high enough pH, 29,30,39 they are relatively unreactive in polymerization, 40 and would add unnecessary complexity to our model. Throughout the remainder of this paper, singly ionized silica monomers are denoted S I . As we will see in section III, the concentrations of S N and S I can be obtained from experimental conditions and the deprotonation equilibrium constant.
We must now specify two additional interactions: those between S I -S I and S I -S N . At first glance, it may be expected that ionized monomers will polymerize in the same way as neutral monomers. However, this is not the case, as the strong repulsion between the negative charges on two S I molecules significantly inhibits the reaction. We assume here, for simplicity, that the polymerization reaction energy is exactly cancelled by the electrostatic repulsion, so that the S I -S I interaction is set to zero. As a result, a system containing only ionized silica and solvent will yield a uniform mixture of mostly monomers and a few oligomers arising from random thermal fluctuations. This is consistent with experimental observations showing that S I monomers formed at high pH by the ionization of monosilicic acid are extremely soluble in water, being responsible for a marked increase in the solubility of amorphous silica at pH>9. 29 As for the S I -S N interaction energy, it should also be attractive, since reactions between neutral and ionized monomers do occur in real systems.
Nevertheless, the enhanced solubility of S I should have an impact on the effective interaction energy, making it less attractive than the S N -S N energy. We have studied the behavior of the model for several values of this parameter (see section IV A) and have obtained qualitative agreement with experimental observations for ε *

Another important difference between neutral and ionized monomers is that interac-
tions involving the latter depend on its orientation -interactions will be stronger or weaker depending on whether a given molecule approaches the negatively charged oxygen or the neutral hydroxyl groups. We take this effect into account by assigning a pointer variable to each lattice site occupied by an S I . This variable "points" to one of the neighbors of that site, to indicate the orientation of the negatively charged oxygen. The remaining neighbors interact with a neutral hydroxyl group.
As a result of distinguishing between the neutral and anionic ends (represented by S -I ) of each S I , we increase the total number of interaction parameters of our model to 6 (the 3 described above plus 3 interactions involving S -I ). Following the same rationale as above, we set the S -I -S -I and S I -S -I interactions to zero. As for the S N -S -I interaction, this corresponds to a reaction of the type: This reaction can be obtained by subtracting reaction (5) from reaction (1). Therefore, the equilibrium constant for reaction (6) is the ratio of the equilibrium constants for reactions (1) and (5). Experimental studies show that the equilibrium constant for the neutral polymerization lies between approximately 20 and 500, 39 while the equilibrium constant for monomer deprotonation is 1.75 × 10 6 (see Ref. 39 and section III A). From these values, the equilibrium constant for reaction (6) is on the order of 10 −4 . This value is very low, so we can safely assume that the reaction does not occur in practice. Therefore, we set the S N -S -I energy to zero.

C. Templated Growth
The final step is to include TPA cations in the model, requiring four additional interaction parameters (TPA-S N , TPA-S I , TPA-S -I and TPA-TPA). We may expect strong repulsion between two TPA cations due to their positive charges. However, this electrostatic repulsion is countered by strong hydrophobic attractions between the propyl groups, which influence the structure of water near TPA cations. 36,41 Recent simulation studies reproduce this hydrophobic effect. 42,43 In tetraalkylammonium ions, electrostatic repulsions and hydrophobic attractions are found to have comparable magnitudes, giving behavior intermediate between those of simple electrolytes and apolar molecules. 36,41,43 In light of this, we have set the interaction energy between TPA cations to zero.
As for the interactions between TPA and silicates, we assume that only the electrostatic attraction between TPA and S -I is important. We assign a relatively strong value of −2 to the reduced TPA-S -I attraction. Van der Waals attractions between TPA and neutral silicates are disregarded for simplicity; these interaction energies are thus set to zero. Indeed, quantum chemistry calculations 37 suggest that the interactions between TPA and neutral silica monomers are of the same order of magnitude as interactions between TPA and water, while the attraction between TPA and anionic silica monomers is much stronger. In section IV C we discuss the effect of using different values for the TPA-S -I attraction. Until now, we have ignored the fact that the TPA cation is much larger than the silica monomer. The ionic radius of a TPA cation is approximately 0.45 nm, 44  to represent this size difference in our model is to include a nearest-neighbor repulsion shell around each TPA molecule. If this repulsion is strong enough, we ensure that all of the first-neighbor sites of the TPA remain vacant, thus effectively increasing its collision diameter. We also include a second-neighbor repulsion between two TPA cations, to eliminate configurations in which two of these large molecules sit too close to each other in a diagonal arrangement on the simple-cubic lattice. The attraction between TPA and S -I and the pointer variable must then be extended to a second-neighbor level. This modification yields a value of √ 3 = 1.732 for the size ratio between TPA and silica monomers in our model, which is a reasonable approximation to the estimate of 1.5 mentioned above.
Most of the calculations presented in this paper deal with solutions containing TPA, which we model using a "large" (or second-neighbor) cation. However, we are also interested in studying the effect of cation size on the properties of the nanoparticles. For this purpose, we have also performed simulations using a "small" (or first-neighbor) cation. This is a reasonable approximation to a tetramethylammonium (TMA) cation, which, at an ionic radius of 0.34 nm, 44 is of comparable size to a silica monomer.
According to the above description, the Hamiltonian of the lattice model is: where superscripts FN and SN denote first-and second-neighbor interactions. The Hamiltonian is written as a sum over all contacts between components on the lattice. Indices 1,  Table I shows the  interaction parameters for the model with only first-neighbor interactions, while Tables II   and III present those for the model with both first-and second-neighbor interactions. For brevity, we refer to these as the first-neighbor and second-neighbor models, respectively. In what follows, we perform most of our calculations on the second-neighbor model.
In summary, the model proposed here assumes that the important parameters controlling silica nanoparticle formation are the initial concentrations, temperature, pH, cation size, and the relative magnitudes of the S I -S N and TAA-S -I attractions. Our model assumes united atom representations on a simple-cubic lattice. In addition, because of the simple cubic lattice the coordination number of silica polymers is increased from 4 to 6; changes in reactivity of silicate species with degree of polymerization 30,31,39 are ignored; and long-range electrostatics are reduced to short-range effective interactions. In spite of these simplifications, we shall see that this model captures the essential physics leading to the spontaneous formation of nanoparticles in the early stages of silicalite-1 synthesis.

III. SIMULATION TECHNIQUES
To probe the behavior and properties of the model described above, we have employed a variety of molecular simulation techniques. Most of these are well documented, and thus we provide only a brief description, focusing on the aspects of particular relevance to the problem at hand. The majority of the results presented in this paper were obtained from Monte Carlo (MC) simulations in the canonical ensemble (N V T ). Parallel tempering and reactive ensemble MC (REMC) were used mainly to complement the N V T results. In what follows, "simulation results" refer to N V T MC unless otherwise specified.

A. Canonical Ensemble Monte Carlo
We have performed N V T MC simulations using a simple cubic lattice with periodic boundary conditions in all three cartesian directions. Two types of MC moves were implemented. The first type of move is a swap, consisting of an exchange in position between two molecules located on different sites. The first site is chosen at random from a list of all 'occupied' sites (a site containing solvent is considered vacant), while the second site is randomly chosen from the global list of N lattice sites. Since a great majority of all sites contain solvent molecules, this procedure allows for significant savings in computer time, by avoiding irrelevant swap attempts. As usual, the move is accepted or rejected based on the Boltzmann factor associated with the configurational change. 45 The second type of move is a rotation of the pointer variable placed on the anionic silica sites. This is implemented by randomly choosing an occupied site and checking to see whether it is occupied by an S I molecule; if it is not, the move is rejected. A new value for the pointer variable is chosen at random from all possible values (6 for the first-neighbor model and 12 for the second-neighbor model). The Boltzmann factor for this trial is calculated and the rotation is accepted or rejected accordingly.
Except where noted, each N V T MC run was started from a random initial configuration, in a cubic simulation box of side L, and allowed to equilibrate for at least 1 million sweeps (usually much more than this). If there are N occ occupied sites on the lattice, then a sweep is defined as N occ attempted moves plus N occ attempted rotations, chosen randomly. The average energy of the system was calculated at the end of each sweep and written to file.
Snapshots of the simulation were produced at evenly spaced intervals during the course of each run. We have also calculated cluster size distributions at different stages of the simulation by implementing the Hoshen-Kopelman cluster counting algorithm. 46 Two silica monomers are considered part of the same cluster if they are connected by a single lattice bond.
A correspondence with experiment was established by setting up the initial configuration based on values calculated from experimental control variables. Typically, the composition of the experimental solutions is given in terms of the mole ratios of the species present.
With this information, the mole fraction of TAA (x TAA ) in the N V T simulation is simply calculated from:  (5). Since the total number of moles is conserved in the equation, the equilibrium relation can be written in terms of the mole fractions of the components: The value of the equilibrium constant (K D ) can be obtained from information on the pK a of the acid involved in the reaction. The pK a of siliceous acids depends on the degree of ionization and may also depend on the degree of polymerization. 30 As explained in section II, we assume that the only ionic silicate present is the singly deprotonated monomer, and thus we adopt the literature value of 9.5 for the pK a of monosilicic acid. 39 In dilute solutions, one can further assume that the concentration of water is essentially a constant at 55.6 mol/l.
Using the pK a for water (14), we obtain a value of 1.75 × 10 6 for K D .
So far, equation (9) has 3 unknowns, so we must introduce 2 other equations. The first is a mass balance on the silicon atom, stating that all silica present in the simulation corresponds to the amount of TEOS added to the solution: The second equation used is a charge balance over all ionic species in the simulation: Solving equations (9-11) simultaneously allows for the calculation of the mole fractions of S N and S I , as well as monitoring the evolution of the pH during the course of the simulation.
Finally, the mole fraction of solvent (x W ) is simply: Note that this procedure assumes that the concentrations of the components are fixed at their equilibrium values. We relax this assumption as discussed below, by performing simulations in the reaction ensemble (see section III C).

B. Parallel Tempering
As we see below, the model proposed in section II is difficult to equilibrate using the conventional Metropolis N V T MC. To improve the efficiency of equilibration, we have employed parallel tempering. 32 This is implemented by performing M independent N V T simulations, running at different temperatures. Each replica is allowed to equilibrate using conventional MC moves (as described in section III A). However, after a preset number of such MC moves a replica exchange is attempted, wherein the complete configurations are exchanged between two replicas running at adjacent temperatures. First, one of the M replicas is selected at random. After this, one of its two adjacent replicas is also chosen randomly (in the case of the highest or lowest temperature, this step is omitted, as there is only one possible choice).
Configurations of the two replicas are then exchanged and the energy difference is computed.
The trial is accepted with a probability given by: As such, the acceptance ratio of the replica exchange depends on the energy difference and on the reciprocal temperature difference between two replicas. Naturally, the closer the spacing between replicas, the higher the acceptance ratio. However, the computer time increases proportionally to M . Therefore, one seeks a compromise involving the minimum number of replicas allowing for high acceptance ratios, to ensure adequate sampling of the desired temperature range.
To achieve this, we have used an iterative procedure to find the ideal value of M and the optimal spacing between replicas. We begin by setting up a simulation containing an initial number of replicas in the temperature range of interest, evenly spaced in 1/k B T . 33 A short parallel tempering simulation is then run, and the resulting average energies and acceptance probabilities are plotted as a function of replica temperature. We then refine the replica temperature grid including more intermediate temperatures where necessary, using the following procedure. The aim of this procedure is to obtain a grid that allows for acceptance ratios of at least 5% throughout the entire temperature range. If exchanges between replicas i and i+1 fall below 5% acceptance, we add more replicas between copies i and i+1, evenly spaced in 1/k B T . To estimate how many replicas to add, we use equation (13) to forecast new replica exchange probabilities. Regarding the energy to insert into equation (13), we use equilibrium internal energies at the new replica temperatures. Since these are not known a priori, we apply linear interpolation to the already-computed internal energies to estimate the new ones. The refined simulation is run and the above procedure is repeated if necessary. In practice, we found that at most three iterations are necessary to obtain the optimal grid. This iterative procedure ensures that temperature regimes with steeper energy variation are assigned a higher number of replicas. We have monitored the movement of replicas during parallel tempering following the procedure described in Ref.

47), and have found sufficient exchange even across regions with significant variations of E
with T .

C. Reaction Ensemble Monte Carlo
As described in section III A, N V T MC simulations are run at a fixed composition, corresponding to the equilibrium composition determined by the silica deprotonation reaction.
An alternative would be to use reaction ensemble Monte Carlo simulations, 34,35 starting from only neutral silica monomers (produced from TEOS hydrolysis), TAA and hydroxyl ions (from the complete dissociation of TAAOH). The main difference relative to N V T MC is that now, equilibration between components related by chemical reactions is sampled (the total number of molecules, however, remains constant). Thus, the concentrations of S N and S I are allowed to fluctuate, which more closely resembles the actual physical system.
To sample these fluctuations in composition, we include an additional type of MC move -a reactive move -which consists of an attempt to transform an S N molecule to an S I (or vice-versa). The probability of accepting a generic reactive move is: where V is the system volume, N i is the number of molecules of component i, ν i is the stoichiometric coefficient of component i (positive for products and negative for reactants), ν is the net change in the number of moles as a result of the reaction (i.e., the sum of ν i over all components), c is the total number of components involved in the reaction, and K eq is the ideal gas equilibrium constant. In order to satisfy detailed balance, the probability of attempting the forward reaction must be equal to that of the reverse reaction.
To ensure that this condition is fulfilled, we proceed as follows. First, a site occupied by an S N molecule is chosen at random and changed to an S I . The energy difference is computed and the trial is accepted with a probability given by: Equation (15) is obtained from equation (14) considering thatν = 0. After the forward reaction trial, a site occupied by an S I molecule is picked randomly and changed to an S N .
The reverse reaction is accepted with probability: Due to the simplifications introduced in our model, we do not explicitly represent hydroxyl ions or water molecules. Thus, the only explicit changes introduced in the composition of the system as a result of the reaction are in the mole fractions of the S N and S I species.
Nevertheless, in the REMC simulations one must keep track of the number of hydroxyl ions that would be present in solution, as this quantity must enter equations (15) and (16) for the acceptance probabilities. The simulation protocol for a REMC run is essentially the same as for N V T MC, except that a sweep is now defined as a swap trial, a rotation trial, a forward reaction trial and a reverse reaction trial per occupied site, chosen with equal frequency.

IV. RESULTS AND DISCUSSION
We begin this section by analyzing the parameter space of the model, identifying the range of conditions that leads to the formation of silica nanoparticles. We then study the stability of these nanoparticles and proceed to compare their properties to experimental observations.

A. Model Parameterization
In this subsection, we examine how changing the interaction parameters introduced above influences the predictions of the model. We begin by studying the effects of the S N -S N energy, considering a system composed of only neutral silica and solvent. This interaction is meant to mimic the driving force for the silica condensation (equation (1)), and so it must be attractive. This attraction will tend to make silica monomers cluster together, while entropic effects will tend to keep monomers randomly distributed in solution. In this case, the system will separate into a phase of pure silica and a phase composed of solvent with some dissolved silica monomers. Just how much monomer remains in solution depends on the ratio between the system temperature and the interaction energy (i.e., on the reduced temperature). In this simplified system, increasing the S N -S N attraction (or equivalently, reducing the temperature) reduces the fraction of monomers in solution, and vice-versa.
We calibrate the value of ε S N S N by comparing simulation results for the simple system  Apart from controlling the charge, the S N -S I attraction also influences nanoparticle structure. In particular, if the magnitude of this attraction is large relative to ε S N S N , nanoparticles show significant sublattice ordering with alternating S N and S I molecules in the solid phase, in a manner similar to antiferromagnets. 48 As the relative S N -S I attraction decreases this effect disappears precipitously. Figure 2 also shows a plot of the percentage of sublattice ordering in the solid phase as a function of ε * S N S I . This was obtained by calculating the average fraction of molecules in the solid connected only to distinct species. Sublattice ordering in the clusters is almost complete for ε * S N S I = −1.2, but decreases rapidly as the energy becomes less attractive. A decrease in the ε S N S I attraction also leads to a higher proportion of S I on the surface of the clusters, rather than inside the core. This is due to the orientational variable assigned to the ionized silica. Because the S N -S -I interaction is set to zero, there is an energetic incentive for S I molecules to be on the surface of a cluster, with the pointers oriented towards the solvent. For ε * S N S I < −0.7, we find no ionized silica inside the core of the particles.
Sublattice-ordered silica clusters seem very unlikely, arising in our simulations as an artifact of the lattice geometry. Therefore, the value chosen for ε * S N S I should produce minimal sublattice ordering. On the other hand, an interaction that is too weak will drive all of the ionized silica to the solution phase. It is known from experiments in colloidal silica 29 that the particles are charged at pH> 7. Silica nanoparticles formed during the clear-solution synthesis of silicalite-1 are also negatively charged. 16 In the remainder of this paper, we have taken ε * S N S I = −0.8, since this value produces little sublattice ordering (less than 1%) and some particle charge. This value also gives the possibility of charge existing inside the nanoparticle core. We note that values in the range of −0.7 to −0.9 also give qualitatively similar results.
Finally, we consider the effect of including TAA cations in the system. Apart from size exclusion, the only relevant interaction involving TAA is the attraction towards S -I . For simplicity, let us consider a solution containing only TAA and ionized silica. Our model predicts that all of the silica remains solubilized with practically no clustering. This behavior arises because we set the S I -S I interaction to zero, and because we use pointer variables for TAA-S -I attractions. Small silica oligomers (dimers and trimers) are occasionally observed in our simulations due to random fluctuations. Such a scenario mirrors experimental studies of TAA-silicate solutions at high pH; these studies report the presence of mostly ionized monomers together with small oligomers such as the cubic octamer. 49 For relatively low values of the TAA-S -I attraction, TAA cations are found in our simulations to be mostly dispersed in the solvent. As this attraction is increased more TAA cations form ion pairs with anionic silica molecules. By setting ε * S -I TAA = −2, most of the TAA and S I ions are paired. Due to the orientational character of the interaction, a given S I molecule is attracted to a single TAA. If the pointer variable were removed, extending the electrostatic attraction to all neighbors of the ionized silica, large sublattice-ordered aggregates of S I and TAA would form in our simulations. These aggregates represent a crystalline phase of silica anions and TAA cations, which is not known to occur at the conditions of interest to this study. Thus, a useful consequence of accounting for the orientational structure of ionic silica is the removal from our simulatioins of S I -TAA clusters.
In summary, a solution containing only ionized silica and TAA cations displays a single phase comprised of isolated or paired ions in solvent, with no significant clustering (see Fig.   3a). This corresponds to an experimental solution where the amount of base is sufficient to ionize all the silica, so that reaction (5) proceeds to completion. As the amount of silica is increased, the concentration of free hydroxyl ions decreases until the equivalence point is reached. Beyond this point the solution has a non-negligible concentration of neutral silica monomers. In situ SAXS and SANS measurements of these solutions show that nanoparticles form once the equivalence point is crossed. 16 Our model reproduces these phenomena. As soon as the mole fraction of S N exceeds the solubility limit, multiple silica aggregates are formed (compare snapshots shown in Figs. 3a and 3b). These aggregates are negatively charged, with the charge residing mainly on the particle surface. Due to the strong attraction between S I and TAA molecules, the nanoparticles are surrounded by several organic cations.
This agrees with the core-shell structure inferred from SAXS and SANS experiments. 16 In the following subsections, we discuss in more detail the nature and stability of these silica nanoparticles, as well as the effects of the model parameters on their properties.

B. Equilibrium vs. Metastability
As described in the previous subsection, when the solution contains a high enough concentration of neutral silica, clusters spontaneously form. It is interesting to examine in detail the evolution of a typical simulation run from a random initial state until cluster sizes reach a plateau. Figure 4 shows a plot of the average energy, average cluster size and average monomer concentration as a function of the number of sweeps in the course of a typical simulation. At early times we observe a steep decrease in the dissolved monomer concentration, because of very rapid formation of small clusters. After about 10,000 sweeps the system essentially runs out of free monomer, containing instead many small silica clusters with an average size of about 60 units.
After this initial stage, the average cluster size increases as the energy decreases. We find that this proceeds by a phenomenon analogous to Ostwald ripening, with dissolution of smaller clusters and growth of larger ones. This is consistent with experimental data showing that silica solubility decreases with increasing particle size. 29 Due to this solubility difference, smaller clusters dissolve, providing monomer for the growth of larger ones. In future work, we plan to study in more detail the dynamics of the growth process. 50 An alternative explanation for the increase in average cluster size is particle aggregation.
We did not observe any cluster aggregation during our simulations. This is consistent with the fact that we did not include global cluster moves in our Monte Carlo scheme. Our clusters move in space only by co-operative sequences of single-particle moves. The high surface charge present in real silica particles at high pH stabilizes them against aggregation through electrostatic repulsion. 29 Therefore, under the conditions studied in this paper, silica nanoparticles likely grow by a mechanism dominated by Ostwald ripening.
Ostwald ripening stops when the difference in solubility between the smallest and largest particles present becomes negligible. After long simulation times, we observe no further Ostwald ripening, no cluster aggregation and no cluster breakup. At this point, we might be tempted to believe that the system has reached equilibrium, characterized by multiple silica clusters in a dilute solution of silica monomers. However, after performing many simulation runs that are thermodynamically identical but statistically independent, we find final states with different average cluster sizes and energies. These properties differ by more than the fluctuations in any given run, showing that each simulation becomes "frozen" in a state that does not necessarily sample from the equilibrium distribution. At this point, two questions may be raised. What is the true equilibrium state of the system? Do the N V T simulations sample well-defined metastable states?
To answer these questions, we made use of parallel tempering Monte Carlo as described in section III B. The freezing of the N V T simulations arises because attractions are strong relative to thermal energy, making it unlikely that N V T simulations can overcome free energy barriers. In contrast, parallel tempering allows the simulation to overcome such barriers, by attempting exchanges between replicas at different temperatures. Figure 5 shows a plot of the average cluster size in a parallel tempering simulation as a function of temperature. A significant finding is that all parallel tempering runs below T * = 0.4 give one single large silica cluster (see also Fig. 3c), whereas all N V T runs below T * = 0.33 give multiple metastable clusters. For T * < 0.4, the mean cluster size from parallel tempering decreases with increasing temperature, reflecting increased silica solubility at elevated temperatures.
In contrast, for T * < 0.33, the mean cluster size from N V T simulations increases with temperature, reflecting the greater likelihood of Ostwald ripening at elevated temperatures. At T * = 0.33, the N V T results jump to the equilibrium curve, suggesting that at this reduced temperature spontaneous fluctuations are sufficient to overcome free energy barriers separating the multi-cluster metastable state from the single-cluster equilibrium state. At T * = 0.4, both N V T MC and parallel tempering predict complete dissolution of all nanoparticles. We note that the transition at T * = 0.33 from a metastable solution of nanoparticles to a single solid is reminiscent of the experimentally observed process of nanoparticles leading to zeolite growth at elevated temperatures. 10,15 Having established that N V T simulations produce metastable nanoparticles, we must now determine whether these metastable states are reproducible. This is crucial for establishing that N V T simulations can be used to model actual nanoparticles, whose metastability has been shown to be quite reproducible from lab to lab. [5][6][7][8][9][10] Towards this end, we have examined average cluster sizes and cluster size distributions (CSDs) from simulations using different system sizes, by performing several independent N V T runs for each lattice size.
After each independent run converged (i.e., the plateau in energy and average cluster size was reached), we calculated the CSD for each realization. Each run required between 2×10 6 and 5×10 6 sweeps to reach convergence, followed by a sampling period of at least 1×10 6 sweeps. A final average over all realizations was then computed. The total number of runs for each system size was determined by keeping the total volume sampled approximately constant from one system size to the next. For example, the calculation for L 3 =60 3 requires 8 times as many runs as that for L 3 =120 3 , since the volume of the latter box is 8 times as large.
A plot of the average cluster size as a function of system size is shown in Fig. 6. For small systems the cluster size increases with L, the linear box size, because a single silica cluster is formed in most of the individual runs. However, for L > 80 multiple clusters are formed in all realizations, and the average cluster size remains practically unchanged. As such, we find that the average cluster size has reached convergence within statistical uncertainty (error bars shown in Fig. 6 represent ± one standard deviation). A more stringent test of this convergence is provided by examining how the CSD varies with system size, shown in can be used to probe the metastable behavior of these nanoparticles.
We have also performed calculations in the reactive ensemble using the protocol described above in Sec. III. The resulting CSD is very similar to that obtained from conventional N V T simulations (data not shown). A characteristic metastable state is also attained in REMC simulations. This lends further credence to the conclusion that this metastable state is representative of the real system, since REMC simulations are a more realistic representation of the experimental situation.

C. Comparison with Experiment
Here we examine the effect of several model parameters on the CSD, and compare trends predicted from our simulations with those found in experiments. 16 We begin by analyzing solutions with different initial concentrations of TPAOH and TEOS, in close analogy with experimental work. 16 Figure 8 shows CSDs computed at room temperature for three different concentrations of TPAOH, keeping the total silica concentration constant. All distributions show the same general shape, but there is a shift in the peak of the distribution to smaller cluster sizes as the TPAOH concentration increases. This means that the size of the most stable silica nanoparticles decreases with increasing amounts of TPAOH. The observation of smaller particles at higher concentrations of TPAOH is in qualitative agreement with experimental SAXS and SANS results, which also show a marked decrease in particle size. 16 Increasing the initial concentration of TPAOH naturally increases the mole fraction of TPA cations in the simulation, but also increases the pH. Therefore, the simulations shown in Fig. 8 are not sufficient to disentangle these separate effects. One can envisage an experiment in which the initial pH is set by a given concentration of TPAOH, but the concentration of the organic cation is increased by addition of TPAX (where X is, e.g., chloride or bromide).
In this way, we test the effect of increasing the TPA mole fraction at constant pH, since Cl − and Br − anions only weakly change the solution pH. In analogy with such an experiment, we have performed an N V T simulation with the mole fraction of ionic silica (x S I ) obtained from the initial concentration of TPAOH, but with the mole fraction of TPA obtained according to: instead of using equation (8). We did not explicitly simulate halide ions in this simulation.
The CSD computed for this situation was compared to that from a simulation at the same pH but lower TPA mole fraction, and to another with the same TPA concentration but higher pH. From the results presented in Fig. 9, we see that increasing the concentration of the organic cations by addition of TPAX has very little effect on the cluster size distribution.
The CSD is almost indistinguishable from the one at equal pH, but very different from the one at higher pH. This suggests that in the range of conditions studied in this paper, it is the pH that controls the stability of silica nanoparticles, and hence determines their size distribution.
The control of particle size by initial pH can be interpreted in terms of nanoparticle structure (see Fig. 3b). This involves a core of mostly neutral silica, together with silica anions located mainly on the particle surface, surrounded by a shell of TPA cations to balance charge. The TPA layer has the effect of protecting the particle surface against further addition of silica monomers. The TPA shell thus contributes to free energy barriers keeping these particles metastable. As the initial pH increases, more of the neutral silica monomers become ionized, so the ratio of S I to S N in the N V T simulation increases. A consequence of this is a significant increase in the particle charge on the surface, and subsequent coverage by TPA cations. Therefore, higher pH means that formation of the protective TPA layer, and hence inhibition of growth, occurs for smaller particles. This effect is similar to a recently proposed stabilization mechanism of silica octamers by TMA cations. 49 TMA cations are seen to adsorb around negatively charged silicates, protecting them against hydrolysis by water molecules. Similarly, TPA cations adsorb around the negatively charged surface of silica nanoparticles, inhibiting further growth by monomer addition.
Our suggestion that nanoparticle size is controlled by surface electrostatics can be exam-ined by performing a sensitivity analysis on our results with respect to the TPA-S I attraction.
In Fig. 10, we show the effect on the CSD from decreasing the value of ε S I TAA by one half.
Reducing this attraction causes a complete collapse of the metastable CSD. Cluster sizes are spread throughout a large size range, and the only distinct peak in the distribution shows up at 3600 units. This peak corresponds to the largest possible cluster given the lattice size, indicating that this N V T simulation has easy access to the global free energy minimum.
These results predict that the attraction between silica anions and organic cations plays a crucial role in stabilizing nanoparticles. Without this stabilizing influence, the simulated particles would grow indefinitely corresponding to formation of a precipitate in experiments.
In Fig. 11, we show CSDs for solutions with different concentrations of silica, keeping the TPAOH concentration constant. The most distinct feature is an increase in the number of clusters present throughout most of the size range (above about 150 silica units). This is in agreement with experiment, which shows that the particle number density increases with increasing TEOS concentration while the particle size remains approximately constant. 16 Our simulations also show a shift to larger sizes. Although this prediction may be in error, it is also possible that this subtle change might be beyond the accuracy of the scattering techniques.
Another experimental observation is the presence of silica nanoparticles in solutions containing other TAA cations. The size of these particles seems to be relatively insensitive to the size of the cation. 16 As explained in section II, our model allows for the study of a smaller cation by eliminating the second-neighbor repulsion around TPA, giving a cation with the same size as a silica monomer. Such a model can be roughly identified with an experimental solution containing TMA ions. A comparison of the results from the first-neighbor (TMA) and second-neighbor (TPA) models is shown in Fig. 12. This comparison shows that the CSD is practically unchanged by varying the cation size. This occurs because under these conditions, the surface density of TAA cations on nanoparticles is still too low for TAA size exclusion to play a role. An additional consequence of this observation is that, contrary to our initial expectations, the simplest model that describes silica nanoparticle formation need not include second-neighbor interactions.

V. CONCLUSIONS
We have developed and applied a lattice model for aqueous solutions of silica and tetraalkylammonium cations at high pH. These solutions represent the early stages of silicalite-1 synthesis, in which silica nanoparticles are formed as precursors to zeolite growth.
Using this model, we observe the spontaneous formation of nanoparticles under conditions at which they are seen experimentally. The nanoparticles are identified as a metastable state, separated by an energy barrier from the true equilibrium state, which is a single large silica cluster. The simulated nanoparticles are negatively charged with the charge residing mainly on the surface, and are surrounded by a layer of tetraalkylammonium cations, in accordance with experimental observations. Starting from a uniform mixture, we observe particle nucleation and subsequent growth by a phenomenon analogous to Ostwald ripening. Nanoparticle growth slows down and eventually stops once a protective layer of organic cations has formed. This layer is stabilized by electrostatic attractions between these cations and the negatively charged silica surface. A similar stabilizing influence of TAA molecules around silicate species has been reported previously. 49 Our model shows excellent qualitative agreement with recent in situ experimental investigations of nanoparticle formation, 16 by reproducing relationships between particle size and experimental variables. In particular, our simulations predict that particle size is mainly controlled by the solution pH. Increasing the pH increases the surface charge, and hence leads to the formation of the TAA layer for smaller particle sizes. At the range of conditions studied in this paper, the concentration and size of the cation have a negligible effect on the size of nanoparticles. Increasing the total silica concentration at constant pH mainly causes an increase in the number of particles formed, but also slightly increases their size.
Increasing temperature causes particle growth due to increased mobility of the TAA layer, until thermal fluctuations are large enough to allow the particles to grow indefinitely towards the equilibrium state.
Due to the coarse-graining of the model, quantitative comparisons with experimental data are difficult to make. Furthermore, the simplifications introduced herein preclude us from predicting the internal structure of the particles. This type of information can best be obtained from continuum atomistic simulations of silica-TAA alkaline solutions. Such simulations may also help clarify the precise role of TAA in nanoparticle stabilization. On the other hand, there is scope for refining the present model; efforts in this direction may yield additional insights into the early stages of zeolite growth.