Monte Carlo Simulations of Silica Oligomerization: A Comparative Study of Parallel Tempering, Umbrella Sampling, and Metadynamics

Abstract

Monte Carlo (MC) simulations of silica oligomerization face significant challenges due to rare event dynamics and high free energy barriers inherent in polymerization reactions. This study introduces a systematic evaluation of advanced sampling algorithms—Parallel Tempering (PT), Umbrella Sampling (US), and Metadynamics (MetaD)—to enhance phase space exploration and accelerate convergence in silica systems. The Beest-Kramer-van Santen (BKS) potential was employed to model interatomic interactions, with explicit parameterization of Si-O bond dissociation energies (ESi-O ≈ 4.5 eV) and angular terms governing tetrahedral coordination. Reaction coordinates such as the degree of polymerization (Qn, where n denotes the number of bridging oxygens) and ring statistics were analyzed to quantify oligomer distributions. PT simulations utilized 32 replicas spanning T = 300–2000 K, achieving exchange probabilities of 15% via optimized temperature spacing. US applied harmonic biases (k = 200 kcal/mol·˚A2) along Qn, while MetaD employed Gaussian deposition (σ = 0.2, ω = 1.2 kcal/mol) every 500 MC steps. Validation against experimental 29Si NMR data revealed PT and MetaD reduced sampling error by 62% compared to conventional Metropolis-Hastings MC. Activation free energies (∆G‡) for trimer formation decreased from 28.3 ± 1.5 kcal/mol (standard MC) to 19.8 ± 0.9 kcal/mol (MetaD), aligning with Arrhenius-derived estimates. Convergence analysis demonstrated PT achieved ergodicity in 106 steps versus 108 for brute-force methods. These results establish that advanced sampling algorithms mitigate kinetic trapping and enable atomistic prediction of silica gelation kinetics under ambient conditions.