Hybrid Monte Carlo-Molecular Dynamics Framework for Silica Polymerization: Bridging Short-Timescale Reactions with Long-Range Structural Ordering

Abstract

Silica polymerization remains a critical process in materials science, yet its multiscale nature challenges conventional modeling approaches. This work presents a hybrid computational framework combining kinetic Monte Carlo (kMC) and reactive molecular dynamics (MD) to simulate silica network formation under aqueous conditions. The kMC module governs stochastic reaction kinetics, including hydrolysis ($\text{Si–O–Si} + \text{H}_2\text{O} \rightleftharpoons 2\text{Si–OH}$) and condensation ($\text{Si–OH} + \text{HO–Si} \rightarrow \text{Si–O–Si} + \text{H}_2\text{O}$), with rates modulated by local proton activity ($a_{\text{H}^+}$). The MD component employs the ReaxFF force field to resolve atomic-scale structural relaxation and stress accumulation ($\sigma_{ij}$) in $\text{SiO}_4$ tetrahedra. A bidirectional coupling protocol synchronizes kMC and MD domains via a weighted Hamiltonian $H = \lambda H_{\text{kMC}} + (1-\lambda)H_{\text{MD}}$, where $\lambda$ adaptively scales with the system’s tetrahedrality parameter ($\langle Q^4 \rangle$). The framework introduces a time-decomposition algorithm to handle disparate timescales, with kMC advancing reaction steps ($\Delta t_{\text{kMC}} \sim 10^{-6}\ \text{s}$) and MD resolving picosecond-scale vibrations ($\Delta t_{\text{MD}} = 0.5\ \text{fs}$). Validation against silica solubility curves ($\log[\text{SiO}_2]_{\text{sat}} = -0.38\cdot\text{pH} + 3.2$) and small-angle neutron scattering (SANS) data reveals consistent oligomer size distributions (radius of gyration $R_g = 1.2–2.8\ \text{nm}$). However, discrepancies emerge in gelation times ($t_g$), with simulations overestimating experiments by 18–22\% at $\text{pH} > 5$, suggesting incomplete treatment of hydrogen-bonding networks. The model identifies metastable $\text{Q}^3$-rich clusters as kinetic traps, delaying percolation thresholds. This work highlights the necessity of explicit solvent modeling and strain-dependent rate corrections in multiscale silica polymerization frameworks.