In this thesis a molecular simulation program is constructed from scratch. It utilizes the Vashishta-Kalia-Rino-Ebbsjø many-body potential model in an eﬃcient parallel implementation. This implementation is validated against physical known properties of silica and water which are known from experimental data. The implementation is then used to investigate a contact problem of a cylinder and half space, both consisting of silica. Eﬀects of surface hydration, and surface roughness are investigated. The results suggest that hydrating silica surfaces will result in lower pressure in the area of contact. In addition the reduced elasticity modulus of a rough surface has a slightly smaller value compared to a smooth surface, according to data gathered from the performed simulations.