Mesoscale model of the synthesis of periodic mesoporous benzene-silica

.


Introduction
Periodic mesoporous organosilicas (PMOs), as most porous compounds, find application in catalysis [1,2], in the adsorption of pollutants [3,4] or gases [5], in controlled drug release [6], among others [7,8].The synthesis of PMOs begins with a process by which the organosilicate species promote the formation of a liquid crystal template from aggregates of surfactants.[9,10] The amphiphilic nature of the latter promotes its self-assembly in an aqueous solution into diverse ordered shapes, such as spherical or prolate micelles, rods, and lyotropic liquid crystals.Surfactants contain hydrophobic tail and hydrophilic head groups, and therefore tend to aggregate in order to keep the tail groups away from the bulk water phase.The alkylammonium surfactants are a class of such compounds, well-known for their application as templates in the synthesis of nanoporous materials.In particular, hexadecyltrimethylammonium bromide (CTAB) is among the most popular surfactants used in the synthesis of PMOs.[11] CTAB comprises a hydrophilic trimethylammonium head and a hydrophobic hexadecane tail.Thus, CTAB forms micelles in aqueous solutions above the critical micelle concentration (CMC), with the positive charge of the ammonium ions at the surface screened by a diffuse layer of Br -ions outside the micelle surface.Experimental studies estimate an average aggregation number of CTAB spherical micelles of around 100 ion pairs.[12][13][14][15][16][17][18] The early stages of the PMO synthesis follow the same mechanism as periodic mesoporous silicas (PMS), with organosilicates as the inorganic moiety instead of silicic acid (TEOS).[11] Nowadays, the most accepted mechanism of PMS formation can be summarized as follows: [19] (i) silicate molecules are added to a solution of CTAB spherical micelles; (ii) silicate molecules replace (at least partially) the bromide ions at the micelle surface; (iii) the silicate moieties at the CTAB micelle surface cause an increase in micellar volume, where flatter surfaces are now more energetically favored; (iv) further silica oligomerization promotes the arrangement of rods into diverse long-range ordered phases such as hexagonal or bicontinuous (a phase formed when several components of the system are intertwined around each other, in contrast with a micellar phase, where one of the phases is arranged in small aggregates [20]); (v) finally, after a silica network is formed around the surfactant template, the latter is removed, usually by calcination [21].This very complex process involves a combination of several simultaneous steps -amphiphilic effects, solvation, thermodynamic equilibrium and chemical reactions among others -which cannot be addressed individually using experimental techniques.Thereby, molecular dynamics (MD) simulations can shed light into this issue, providing a microscopic point of view of how these steps occur and how they can be controlled to produce tailor-made PMSs for specific applications.[22][23][24][25] The first ordered mesoporous organosilica hybrid material with a crystal-like wall structure used benzene as the organic linker (the organic moiety that connects the inorganic groups) between two silicate moieties, [26] and the synthesis procedure was reported to be as simple as adding benzene-bridged organosilica monomers to a 3.3 wt% aqueous solution of alkyltrimethylammonium surfactants, followed by removal of the surfactant through solvent extraction.The powder X-ray diffraction patterns of the as-made material containing surfactant and the surfactant-free (periodic mesoporous) benzenesilica (BZS) revealed that both were arranged in two-dimensional hexagonal lattices, with lattice constants differing by less than 2%.Transmission electron microscopy (TEM) images and electron diffraction patterns indicated that the resulting material was composed of one-dimensional hollow rods, with diameters of ca.38 Å, arranged in a hexagonal manner.However, the role of the organosilicates in the formation of PMOs remains poorly characterized.In particular, it is not clear whether the organosilicates simply condense around preformed surfactant mesophases, or if they are the driving force for the transition from individual micelles to the hexagonal rod array.In the PMS synthesis, the latter has been shown to be the case.[23] Moreover, from the CTAB phase diagram, it is known that, at room temperature and a concentration of 3.3 wt%, the surfactants should form mostly spherical aggregates.[25] This suggests that BZS could play an active role in the PMO synthesis from the very initial stages, in the so-called cooperative templating mechanism.[20] Futamura et al. have successfully performed density-functional theory (DFT)-based calculations to predict the most stable conformations of benzenesilicates, [27] and to calibrate the parameters of an all-atom (AA) model of organosilicate precursors.Subsequent AA MD simulations (with the smaller decyltrimethylammonium as template) yielded surfactantorganosilicate spherical micelles.[24] In previous simulations of PMS synthesis, the bromide counterions were exchanged by negatively-charged silicates, which adhered to (and slightly penetrated) the micelle surface, promoting micelle growth and further aggregation.[22,23] The AA simulations of surfactant-organosilicate micelles revealed that the role of the negativelycharged organosilicas could be very similar in the synthesis of PMOs.[24] However, the system size and time scale limitations inherent to AA simulations only allow the simulation of a rather small number of surfactants, which hampered the study of mesophases such as the lyotropic hexagonal phase observed experimentally.[26,28] For example, a recent AA md study was able to model a pre-shaped hexagonal phase of PMO, but not the dynamics that lead to this phase, due to limitations of simulation time and system size.[29] Much more simplified lattice models have been applied to further understand phase diagrams of hybrid organic-inorganic systems, [30][31][32] but here the problem is to establish a direct connection between the results of the model and realistic experimental systems.The use of a bottom-up multi-scale modelling approach, based on developing a mesoscale coarse-grained (CG) description of the system, can overcome both of the above limitations, allowing the use of larger systems and much longer simulation times than in AA models, but maintaining the necessary degree of chemical specificity to enable a direct link with experimental systems.[33] In the CG scenario, the system is not modeled as individual atoms but instead as interaction centers (called beads) representing groups of atoms.This considerably reduces the number of interactions and allows the simulation of larger systems for longer times.The MARTINI CG model [34] considered in our work establishes a four-to-one heavy-atom mapping.For instance, four water molecules are modeled as a single interaction center, representing a decrease in the number of distinct particles to be treated by a factor of 12. Effectively, using the CG method, one can model systems with a number of particles about one order of magnitude higher than in AA simulations.
In this work, we developed a CG model based on the MARTINI matrix of interactions to mimic the initial stages of PMO templated synthesis that reproduces and significantly extends the results of previous AA simulations.[24] Aqueous solutions of CTAB and BZS under different degrees of deprotonation were considered to represent the pH conditions used in the experimental synthesis of the phenylene-PMO.[11] In particular, the calibration of our CG model was carried out by comparing the radial density profiles of the micelles obtained by the AA and CG levels of description, as done in previous studies.[23] Different degrees of BZS deprotonation (neutral, singly-charged or doubly-charged) were considered.CG models based on the MARTINI force field have demonstrated the ability to reproduce surfactant phase behavior.[34] Pérez-Sánchez et al. applied this method to model the initial stages of the PMS synthesis, [23] which involved the introduction of two new CG beads to model the charged and neutral tetraethyl orthosilicate monomers (TEOS).After a careful calibration of our CG model for BZSs, as described in more detail below, MD simulations containing a few thousand CTABs and hundreds of thousands of water molecules were considered (corresponding to a concentration of 13 wt%).The evolution of the system towards the mesophases observed experimentally, i.e., the sphere-to-rod transition and the subsequent hexagonal rod arrangement, was analyzed.In the remainder of this paper, the computational details and CG model are disclosed in section 2, the results and discussion are presented in section 3, and finally section 4 summarizes the conclusions.

Simulation Details
All simulations were carried out with GROMACS 5 [35] using a leapfrog algorithm [36] to integrate the equations of motion.The time steps for the AA and CG simulations were 2 fs and 10 fs, respectively.Orthorhombic simulation boxes with periodic boundary conditions in all cartesian directions were used.This box shape does not affect the resulting phases of the system, [19,25] and was chosen to speed up the formation of ordered mesophases.In both AA and CG simulations, the temperature was maintained at 298 K (the temperature at which the benzenesilicate precursor and the surfactant solution are mixed experimentally [26,28]) using the velocity-rescaling thermostat, [37] while the pressure was fixed at 1 bar with the a i ll − a a ba tat.[38] The energy contributions to the potential energy function were the bond stretching, angle bending, and torsion for bonded interactions, whereas the L a d− L ) a d Coulombic terms were considered in the non-bonded interactions.The AA models were taken from our previous works, [22,27] and detailed tables of parameters are provided in those publications.The rigid SPC/E potential [39] was used to represent the AA water molecules.Bond lengths were constrained by the LINCS algorithm.[40] In the AA model, short-range dispersion interactions were handled using a twin-range cutoff scheme, with interactions between the inner and outer radii of 1.0 and 1.2 nm, respectively.A potential cutoff with the force-switch modifier function was employed in the CG approach, using a cutoff radius of 1.2 nm where the energy decays smoothly to zero from 1.0 to 1.2 nm.Longrange Coulombic electrostatic interactions were evaluated with the Particle Mesh Ewald (PME) method [41] employing a cutoff radius of 1.2 nm.

Coarse-grained Mapping and Model Development
To model aqueous PMO synthesis solutions, we need to describe several components in our CG model: water, surfactants, counterions (bromide and sodium), and organosilicates.The CG mapping for the CTA + surfactant was validated and successfully used in previous studies.[23,25,[42][43][44] As shown in Figure 1, the CTA + model consists of four apolar C 1 beads (see reference [34] for bead nomenclature) to describe the alkyl-chain tail and the charged Q 0 CG bead represents the headgroup.Bromide and sodium ions are represented by charged beads with acceptor and donor hydrogen bonding capabilities, respectively (i.e., Q a and Q d ), both of them implicitly solvated by six water molecules as per the MARTINI specification.[34] Each CG water bead, P 4 , reproduces the behavior of four real water molecules.In all the CG simulations, 10% of slightly bigger antifreeze water beads, BP 4 , were included to avoid the tendency of the P 4 water to unrealistically pack into a lattice and freeze.[34] We note that the antifreeze particle only interacts differently with other water molecules, while interactions with different species are kept identical to prevent any disruption to the system.The BZS molecule is composed of an aromatic ring with two substituted silicate groups in para position (Figure 1).The standard MARTINI model makes use of three smaller SC beads to represent aromatic rings.[34] However, for para-substituted rings, this mapping would modify the required symmetry of the molecule, as illustrated in Figure S1.Hence, a different approach was adopted here, whereby each BZS molecule was composed of four CG beads in a diamondlike arrangement: two beads centered on the silicates at the end of the molecule and two beads to represent the benzene organic ring.This choice is justified by the need to maintain the original structure and symmetry of BZS.Because we changed the way we described the aromatic ring, relatively to the original MARTINI model, we adjusted the interaction potential of those beads to best reproduce the AA micelle density profiles, as discussed below.The silicate moieties are known to protonate and deprotonate in response to changes in the pH of the medium.[45] Because each BZS molecule contains two silicate groups, and assuming, as we have done in the past, [25,46] that each silicate can only lose one proton (i.e.there is a maximum of one charge per silicon atom), we need to consider three distinct protonation states: IBI represents a BZS molecule where each of the two silicate groups is deprotonated, hence negatively charged; NBN represents a molecule with two neutral silicate groups; IBN represents a molecule where only one of the silicate groups is negatively charged.Experimentally, the synthesis of PMOs is commonly carried out under relatively high pH values, corresponding to most of the BZS being deprotonated.[24,28] In fact, under the representative experimental pH of 11, and assuming the acidic constant of BZS is the same as that of the monomeric uncharged silicic acid, [45] the IBI:IBN:NBN ratio is around 1:0.03:0.0009.Nevertheless, the three versions of BZS were modelled here for generality and to enable simulations to be run at variable pH conditions.The P SN and Q SI CG beads previously developed for the protonated and deprotonated tetraethyl orthosilicates (TEOS), [23] respectively, were used to describe the silicate moieties of the BZS.Thus, NBN comprises two P SN CG beads whilst IBI is made up of two negatively charged Q SI beads.The singly deprotonated species, IBN, includes one of each silica bead type.
Simulations with individual pre-formed spherical micelles built with Packmol [47] were used to calibrate and validate the CG model.Afterwards, simulations starting from randomly arranged initial surfactant configurations were developed to study the CTAB self-assembly in the presence of BZS.The pre-formed micelle simulations were performed at both AA and CG levels while the randomly arranged systems were studied only with the CG model to probe whether the proper micelle distribution was attained.Three pre-formed AA and CG micelles were built, containing: fully protonated NBN, partially deprotonated IBN and fully deprotonated IBI.The number of compounds of each type considered in the AA and CG simulations of one micelle is reported in Table S1.The number of surfactants was fixed at 100 for these simulations, to match the experimentally measured aggregation number of CTAB (i.e.so that only one micelle, on average, is expected to form).The number of water molecules was fixed at 50000 for all charge states of the BZS, corresponding to about 4 wt% concentration of surfactant, therefore well above the CMC threshold.All systems were simulated for 1 s to ensure that equilibrium was attained.The number of BZS included was fixed at 50 in the CG simulations and the number of Br -was selected to balance the net charge as shown in Table S1.Bromide ions were added in the cases where the BZS charge did not completely balance that of the CTA + .In the IBI, IBN and NBN systems, the pre-formed micelle density profiles were obtained with an in-house cluster counting code [48] ba d t H −K p l a cluster counting algorithm.[49] The criterion adopted was that two individual CTA + tails belong to the same cluster when the distance between the last two beads in the chains was lower than 1.2 nm.
The first attempt to parameterize the organic ring was done using the original SC 3 parameters from the MARTINI model.[34] However, this choice led the BZS molecules to penetrate too deeply into the micelle instead of lying at the surface, as predicted by the AA calculations.Indeed, the radial position of the BZS critically depends on its relative interaction strength with water (P 4 ) and the tail atoms (C 1 ); the overall interaction of BZS with the former should be slightly more attractive than with the latter.In fact, the relatively small SN 0 half polar/half apolar bead with no donor/acceptor hydrogen-bond capabilities captures this relative difference quite well.Furthermore, the SN 0 GC bead also takes into account the close proximity of the aromatic beads to the terminal silicate groups, which confer an additional polar character to the molecule.[27] The best agreement was indeed found in the three systems when the SN 0 bead was chosen, as illustrated in the density profiles shown in Figure 2. The parameters for the BZS CG beads used in this work are shown in Table 1.
The CTAB/NBN/W, CTAB/IBN/W and CTA + /IBI/W CG systems were additionally run by placing all the molecules in random positions at the beginning of the simulation, to verify that the proper micelle size was found after equilibrium was achieved (Figure S2).The density profiles closely resembled those obtained using the pre-formed micelles (see Figure S3 and Table S2).All systems aggregated into a single spherical micelle (comprising all 100 surfactants) after approximately 70, 150 and 190 ns, for IBI, IBN and NBN, respectively.In the case of IBI and IBN, all BZS were located at the micelle surface, swapping with the Br -counter anions, while about 80% of the NBN remained solvated.This can be mainly attributed to the lack of Coulomb attraction of NBN towards the surfactant head groups.
Table 1.Lennard-Jones well-depth parameters, ε, in kJ/mol, for the CG benzene ring particle (SN 0 ), and for the ionized (Q SI ) and neutral (P SN ) CG silicate moieties of the BZS.The capital letters between parentheses denote the interaction levels according to the MARTINI notation.The value of the LJ site diameter is  = 0.47 nm for the interaction between all types of particles, except between S particles, for which  = 0.43 nm a d ε i al d t 75% f t di play d valu .The present CG model exhibited a very good agreement with the AA density profiles, reflecting a similar picture in both descriptions.[44] In fact, all profiles reflected the same micelle composition: a hydrophobic core containing exclusively tail atoms, practically isolated from the solvent by a layer of head groups, with the layer of bromide counterions located slightly farther from the center of the micelle.The CG density of bulk water is underestimated by around 10% due to the approximate nature of the coarse-grained water model.[34,50] Gaussian functions were fitted to the normally distributed species (i.e., all except the tails and waters) to provide a more detailed insight.Table S2 shows the relative percentage difference between the abscissas of the CG and the AA peaks.The values in Table S2 tell us that, in general, our CG model produces micelles around the same size as the AA ones.These results suggest that the beads introduced in reference [23] for the silicates are transferable to Journal Pre-proof other models, namely to organosilicates.In addition, the density profile for IBI resembles that of charged silica monomers, as found in reference [23], in the sense that the IBIs lie closer to the center of the micelle than the heads do, contrarily to the case of silica dimers.This means that IBIs, composed of two charged silica moieties, behave like two silica monomers rather than a dimer from a structural point of view.This is possibly due to the larger distance between the negatively charged particles (compared to the Si-Si distance in a silica dimer) and because of the somewhat hydrophobic nature of the organic linker.

Results and Discussion of PMO Self-assembly Modelling
The CG model developed in this work opens the door to the exploration of the initial stages of PMO synthesis, emulating the experimental steps described in the literature, from micelles to the final hexagonal arrangement of rods through sphere-to-rod transitions.Previous experimental works on PMO synthesis [28] p t d t at t Z wa add d " l wly" to a CTAB system formed by spherical micelles.Since the experimental PMO synthesis is done under high pH conditions, unless otherwise stated, all the simulations were carried out with the IBI BZS.Bearing this in mind, CTAB + IBI + Water systems at a surfactant concentration of 13 wt% were selected.Although this concentration is higher than the one commonly used in the experimental PMO synthesis (around 3 wt% [4]), it is still within the CTAB micellar region [25] and decreases the demand of computer resources associated with larger simulations boxes that would be required for modeling more dilute solutions.
To explore the role of IBI in the initial stages of the PMO synthesis, three types of simulations were performed: (R) all the IBIs and CTAB randomly placed in the simulation box, (P) IBIs added in one batch to a previously relaxed solution of CTAB spherical micelles, and (D) IBI add d "d pwi " t a p -equilibrated CTAB solution.Preliminary tests showed that an IBI/CTAB ratio of 1/2 is enough to obtain hexagonally packed rods, while ratios higher than this l ad t t a fi al tat f t y t wit t " xt a" I I lvated (see Figure S4 in the SI).Furthermore, and as shown in previous computational studies, [25] the formation of parallel rods requires a certain minimum number of surfactants.Thus, systems of three different sizes (1000, 2000 and 3000 CTAB) were prepared, while keeping the BZS/CTAB ratio fixed at 1/2.Table 2 shows the number of particles present in each simulation.Systems were labelled using the number of CTABs (in thousands) followed by a letter, indicating random initial positions for all components (R), IBIs added to a pre-equilibrated solution of CTAB micelles (P), or slow dropwise addition of IBIs to pre-equilibrated CTAB micelles (D).In Figure 3, the final MD simulation snapshots with 1000 CTABs are shown.The simulation time was 3 µs for runs 1R and 1P, and 13.5 µs for 1D.In the 1R simulation, the final configuration includes an aggregate of rods and spherical micelles in a square lattice arrangement and two cross-shaped rods, while simulation 1P successfully formed hexagonally arranged rods after 1.4 µs.A visual inspection of the rod diameter displayed a value around 38-41 Å, in good agreement with the experimental value of 38 Å. [26] Some relevant intermediate states of the system during the formation of rods and the hexagonal packing in simulation 1P are shown and described in Figure 4. To ensure that the system attained the proper thermodynamic equilibrium, the IBI molecules (and their counterions) were removed from the final 1P system.Then, the system was equilibrated again and, after 1 µs, the rod-like structures had reverted to spherical micelles (see Figure 5).The 1R and 1P results demonstrate that the PMO hexagonal arrangement can only be attained when the BZSs are added into a CTAB micellar solution rather than simply mixing all the components in a disordered fashion.It is clear that the role of the BZS is important from the very initial stages of the synthesis.The 1P simulation additionally allowed us to conclude that, at the temperature and concentration considered, the CTAB micelles are unable to transition from spheres to rods without the aid of IBIs, regardless of the size of the system and simulation time.Although it is beyond the scope of this work, it is worth mentioning that an estimate of the free energy might allow us to distinguish between a kinetically trapped structure and a thermodynamically most stable one.

Initial
+20 ns +4 ns +85 ns +300 ns +1 µs  The 1D simulation mimics the experimental steps commonly followed in the PMO synthesis.Figure 3 shows the final simulation snapshot of the 1D system depicting the initial hexagonal rod formation of PMO synthesis.The BZS dropwise addition procedure involved performing 8 consecutive simulations, each adding 50 IBIs (and 100 Na + counterions, for a total of 400 IBIs and 800 Na + ).After each addition step, the simulations were run until an equilibrium state was reached.The final configurations after each IBI addition step and further equilibration are shown in Figure 6.The equilibrations with 50, 100 and 150 IBI lasted 1 µs each.Two main differences can be observed between the final state of these three simulations: first, the amount of solvated Br -increases as they are replaced by IBI moieties on the surfaces of aggregates; second, the number of aggregates steadily decreases as they become more prolate-shaped ones.Along the equilibration with 200 IBI, an infinite rod becomes fully formed after 0.6 µs.For this reason, the total simulation time in this step was increased to 6 µs, after which two long rods were found.The 1:5 IBI:CTAB ratio is, therefore, high enough to allow for the formation of long rods, but not for their hexagonal arrangement.After the addition of another 50 IBI and equilibration for 1 µs, no relevant changes were observed.However, as soon as the system contained 300 IBIs, after 1.5 µs, it adopted a configuration with all CTAB + IBI arranged in rods, albeit separated by a considerable distance.Afterwards, the system with 350 IBIs displayed close contact between rods and finally the hexagonal arrangement was achieved with 400 IBIs, as illustrated in Figure 6.Note that the final hexagonal arrangement obtained in 1D is quite similar to the result obtained in the 1P simulation.Thus, from a CTAB micellar solution, the same thermodynamic equilibrium can be attained regardless of the way in which the BZS is added.For this reason, the simulations with 2000 and 3000 CTAB were only carried out starting from random positions (R) and from preformed micelles (P).

Add IBI/Na +
Remove IBI/Na + Journal Pre-proof Increasing the system size in the 2R/P and 3R/P systems (respectively twice and three times as many CTABs and IBIs as included in 1R and 1P) yielded similar results.The addition of IBI to preformed CTAB micelles resulted in hexagonally arranged rods (four rods, in the case of simulation 2P, as shown in Figure S5), while randomly arranged CTABs and produced a disordered arrangement.
An important synthesis stage observed in the above simulations was the sphere-to-rod transition induced after the addition of IBI, resembling the same stage observed in PMS synthesis.[23] Figure 7 shows some selected snapshots of the evolution of the 1P system at different simulation times, displaying the fusion of two spherical micelles of around 100 surfactants each (Figure 7a). Figure 7b shows the first contact between the two micelles with the IBI acting as a bridge between them.In Figures 7c-e, a wider IBI bridge linking the two micelles was formed.Figure 7f illustrates how some surfactant tails (green) of the two micelles come into contact, forcing the system to relax towards a more stable rod-like shape as shown in Figures 7g-h.This role of t Z a ti g a "i i b idg " b tw ig b aggregates to promote sphere-to-rod transitions, resembles the role of silicates in the MD simulations by Pérez-Sánchez et al. [23] devoted to the formation of PMS materials.The pH effect in the PMO hexagonal phase formation was also analyzed for the 1P system, starting from pre-equilibrated CTAB micelles with 250 IBIs and 250 IBNs added in random positions, corresponding to a pH of around 9.5.[45] The resulting configuration was equilibrated for 3 µs and is shown in Figure S6.After 2.2 µs, the system became essentially composed of three hexagonally packed rods, very similar to the those found in simulations 1P and 1D.Therefore, the decrease in pH simply increased the time required to form hexagonally arranged rods from 1.4 µs to 2.2 µs.Thus, the role of a higher pH is, as expected, to promote the deprotonation of the BZSs which connect neighboring rods to form the hexagonal arrangement and, consequently, to accelerate the PMO synthesis process.Figure S7 shows snapshots of one of the first micelle fusions that occurred during the simulation.At this stage, the number of adsorbed IBIs on the micelle surface was noticeably larger than that of IBNs, so that the first contact between the micelles was mediated by surface IBIs.However, as the number of adsorbed IBIs became similar to that of adsorbed IBNs, the bridge between the micelles became composed of both types of BZS, in approximately equal shares, so that both of them contributed to the fusion process.Note that, for the fusion of the two micelles in question, the time scale is roughly the same as the one shown in Figure 7 for the CTAB + IBI system.

Conclusions
In summary, we have developed a CG model to computationally reproduce the early stages of the formation of periodic mesoporous organosilicas from benzenesilicate (BZS) precursor species.The CG mapping of the BZS molecules consisted of four beads: two beads to model the ring, and two silica beads taken directly from previous works on silicic acid.[23] The model was calibrated by reproducing the correct radial density profile for 100-surfactant (CTAB) spherical micelles, when compared to all-atom calculations.In fact, a bromide-to-BZS exchange on the micelle surface was observed after BZN addition to a CTAB micellar solution.The CG micelle density profiles match the AA ones for neutral, singly-and doubly-charged benzene-silicates, supporting the transferability of the special silica beads developed in our previous works.
Three system sizes were simulated, with 1000, 2000 and 3000 surfactants and doubly deprotonated BZS (denoted IBI) with an IBI:CTAB ratio of 1:2.Three types of simulation were compared: randomly adding all the components (IBI and CTAB) at the same time, randomly adding all IBIs to a pre-equilibrated solution of CTAB micelles, and adding IBIs in small steps (dropwise) to a solution of CTAB micelles.Overall, the systems successfully yielded hexagonally packed rods when the IBIs were added to a pre-equilibrated solution of CTAB micelles.In contrast, when both components were randomly added at the same time, amorphous configurations were obtained.This may suggest a phenomenon of kinetic trapping, explaining why in experimental synthesis the silica source is always added to a pre-existing micellar solution.
Our simulation of the dropwise addition of IBIs to a CTAB solution resulted in a configuration similar to the one found from the addition of all IBIs at the same time, suggesting that the hexagonal arrangement of rods is always attained as long as the surfactant solution is pre-equilibrated.In addition, removal of the BZS from a PMO system, leaving only the CTAB, causes the system to revert to spheres demonstrating that thermodynamic equilibrium was reached.In this regard, however, any comparison to experiments is merely qualitative -this is because each of our simulations takes place at a fixed silicate speciation, whereas silica polymerization reactions are likely to take place in experiments quite soon after mixing the two solutions.It is possible that, in PMO solutions, the reaction takes place on a time scale that is comparable to the silica/surfactant self-assembly, hence a dropwise addition is needed to avoid the system becoming trapped in less ordered states.Further work, possibly using a reactive model of silica, [51] is needed to fully address this question.
Even though, at the pH values of 11 typically employed in PMO synthesis, around 97% of the BZS are expected to be doubly deprotonated, we simulated a lower pH, around 9.5, by including a fraction of singly deprotonated BZS, and found that the final configuration of the system is nearly the same as when only doubly deprotonated BZSs are added.In both situations, micelle fusion was found to be mediated mainly by IBIs and to occur on approximately the same time scale, although the hexagonal arrangement of rods took roughly 50% longer to occur with lower pH.When both IBIs and IBNs are present, the IBIs assume the main role in bridging two neighbor micelles and, in later stages of the fusion, some IBNs participate in the sphere-to-rod process.
Contrary to the PMS synthesis, our MD results suggest that silica oligomerization seems unnecessary in the initial stages of the hexagonal arrangement of PMO formation.This is most likely because the organic linker is already connecting two silicate groups, therefore the BZS molecules are able to act as bridges between aggregates even without polymerization reactions.According to our simulations, rods are only formed if the number of added BZS is at least around a fifth of the number of CTABs.Furthermore, the hexagonal arrangement is only achieved above the IBI:CTAB ratio of 2:5.The formation of rods occurs around 20 ns of simulation time after the addition of BZS, and their hexagonal arrangement ensues around 280 ns later.Encouragingly, our estimated rod diameters of 38-41 Å match the experimental value of 38 Å.The detailed mechanism via which rods are formed from spheres and arranged into hexagonal mesophases was found to be analogous to that of PMS, with the aforementioned differences arising from the inherent linkages between monomeric organosilicate precursors.This shows that, similarly to the formation of periodic mesoporous silica, the benzenesilica plays an active and cooperative role in the synthesis of periodic mesoporous organosilica.The special BZS geometry and the IBI density charge resembles the effect of the silica TEOS oligomerization in PMS synthesis.Further research on this matter may include using our modelling approach to describe systems that contain other linkers, such as ethylene or ethane.

Figure 1 .
Figure1.Scheme of the coarse-grained mapping considered in this work.The hexadecyl CTAB alkylchain tail is modeled by four apolar C 1 beads (green) and the headgroups are mapped by a Q 0 bead (purple).[23]The monoatomic ions Br -and Na + are represented by a Q a (teal) and a Q d (orange) particle, respectively (note that these two particles implicitly include six water molecules).The CG water bead in MARTINI, P 4 (blue), includes four atomistic water molecules.Finally, the benzene rings are represented by two SN 0 (grey) particles, and each silica moiety of the benzene-silicates is modeled by a Q SI (deprotonated) or P SN (protonated) particle (red).

Figure 2 .
Figure 2. Radial density profiles of the pre-formed spherical micelle systems used to calibrate the BZS model for IBI (top), IBN (middle), and NBN (bottom).Solid and dashed lines represent the AA and CG data, respectively.The green and magenta lines correspond to the CTA + tails and heads, respectively, blue to water, cyan to Br -counterions, and red and black to the silicate and the organic rings of BZS, respectively.

Figure 3 .
Figure 3. Final state of the simulations with 1000 CTABs: 1R (left), 1P (center) and 1D (right).The color code is as follows: CTAB alkyl-chains and headgroups are depicted in green and purple, respectively, whereas the silica moieties of IBIs are shown in red and their benzene rings in grey.Water molecules and counterions were removed for clarity.

Figure 4 .
Figure 4. Snapshots of a 3 µs-long simulation of an initially pre-equilibrated 1000-CTAB solution of spherical micelles, plus 500 IBI initially placed in random positions (top left).The simulation time difference with respect to the previous snapshot is indicated below each image.Most of the IBI quickly adhere to the micelle surfaces, replacing the Br -, and long rods are formed (top right, Br -are solvated and have been omitted for clarity).The long rods slowly incorporate the remaining spherical micelles and ultimately arrange themselves hexagonally (bottom right).The system retained this state for the remaining 1.6 µs of simulation.Water molecules are hidden for clarity and the color code is the same as in Figure 3.

Figure 5 .
Figure 5. On the left, a 1000-CTAB (green) solution of spherical micelles.At the center, the state of the system after adding 500 IBI (red) and 1000 Na + and equilibrating for 3 µs.On the right, the configuration obtained by removing the IBI/Na + again and equilibrating for 1 µs.Water molecules are hidden for clarity and the color code is the same as in Figure 3.

Figure 6 .
Figure 6.Snapshots of a simulation (lasting 13.5 µs total) f t xp i tal " l w" increments of 50 IBI) addition of 400 IBI molecules to a pre-equilibrated 1000-CTAB solution of spherical micelles.The top left image shows the pre-equilibrated CTAB micelles with 50 IBI placed in random positions.Each subsequent image corresponds to the state of the system after being initially equilibrated (top center), after the inclusion of an additional 50 IBI for a total of 100 IBI and equilibrating again (top right), etc., as

Figure 7 .
Figure 7. Snapshots for the CTAB + IBI + H 2 O system with 13 wt% surfactant concentration, illustrating in detail the sphere-to-rod transition.The simulation time in nanoseconds is shown at the bottom of each frame.Water molecules are hidden for clarity and the color code is the same as in Figure 3.

Table 2 .
Number of particles present in each simulation.The letters R, P and D correspond to random, pre-equilibrated and dropwise simulations, respectively.