Polymerization of 216 silicic acid molecules, H4SiO4, has been studied using the molecular dynamics computer simulation technique. Multibody potentials which reproduce bulk and surface structures of silica and silicate glasses, as well as molecular configurations containing Si, O, and H ions, were used in the simulations. Results of the simulations are consistent with experimental data of sol-gel systems. Chains form at the early stages of polymerization, followed by ring formation, consistent with interpretations of NMR data and semiempirical molecular orbital quantum calculations. The activation energy for formation of branching Q(n) species is 12 kcal/mol, consistent with the experimental data of 12 kcal/mol for gelation. The relative time evolution of the various Q(n) species is consistent with experimental NMR data.