We have used Monte Carlo simulations to study the formation of the MCM-41 mesoporous silica material, with a new lattice model featuring explicit representations of both silicic acid condensation and surfactant self-assembly. Inspired by experimental syntheses, we have adopted the following two-step "synthesis" during our simulations: (i) high pH and low temperature allowing the initial onset of mesostructures with long-range order; (ii) lower pH and higher temperature promoting irreversible silica condensation. During step (i), the precursor solution was found to spontaneously separate into a surfactant silicate-rich phase in equilibrium with a solvent-rich phase. Lamellar and hexagonal ordering emerged for the surfactant-silicate-rich mesosphases under different synthesis conditions, consistent with experimental observations. Under conditions where silica polymerization can be neglected, our simulations were found to transform reversibly between hexagonal and lamellar phases by changing temperature. During step (ii), silica polymerization was simulated at lower pH using reaction ensemble Monte Carlo to treat the pH dependence of silica deprotonation equilibria. Monte Carlo simulations produced silica-surfactant mesostructures with hexagonal arrays of pores and amorphous silica walls, exhibiting Q(n) distributions in reasonable agreement with Si-29 NMR experiments on MCM-41. Compared with bulk amorphous silica, the wall domains of these simulated MCM-41 materials are characterized by even less order, larger fractions of 3- and 4-membered rings, and wider ring-size distributions.