A novel Hybrid Monte Carlo simulation method is used in conjunction with an accurate pairwise additive interatomic potential (Tsuneyuki et al., 1988, Phys. Rev. Lett. 61, 869) to study the phase diagram of silica. In the Hybrid Monte Carlo formalism discussed in this work trial moves are proposed by resorting to molecular dynamics techniques. However, since hybrid methods do not suffer from the numerical instabilities associated with finite step sizes, we are able to use effectively time steps that are several times larger than those required for conventional molecular dynamics simulations. More importantly, such methods provide a means of conducting simulations of complex systems within a Monte Carlo framework, thereby facilitating the study of isothermal or isobaric ensembles and that of open systems and phase equilibria. To illustrate the applicability and the usefulness of the methods we use them to reproduce the various thermal- and pressure-induced solid-solid phase transitions undergone by silica over wide ranges of pressure and temperature. This is of particular interest for mineralogical research at extreme conditions such as those encountered in the Earth’s interior.