Molecular dynamics computer simulations were used to study surfaces of pure silica glass. The potentials used here were those previously established to model bulk silica and have been extended to study surface relaxation in a perfect vacuum. A large number of surfaces were made using different starting configurations; system sizes, and cooling procedures. Following “fracture”, many broken bonds rearranged in response to the changes in the net forces in the surface region. After this reconstruction, the simulations showed the expected general features observed experimentally, such as a prevalence of oxygen atoms at the outermost surface, non-bridging oxygens, and strained siloxane bonds. Three fold silicons (similar to e’ centers) were initially present in the “fractured” surfaces but most often were incorporated into the network tetrahedrally after reconstruction. Other defects produced during the reconstruction were five coordinated silicons and more importantly, edge sharing tetrahedra, forming the strained siloxane bonds. Bond angles and bond lengths for each defect were determined, showing good agreement with previously published results as well as providing new information. Finally, estimations for silanol concentrations were made which compare well with experimentally determined coverages. The computer simulation technique used here adequately reproduces many of the structural and dynamic characteristics of silica glass surfaces.