Shadow Hamiltonian in classical NVE molecular dynamics simulations: A path to long time stability

The shadow energy, E-s, is the conserved quantity in microcanonical ensemble (NVE) molecular dynamics simulations carried out with the position Verlet central-difference algorithm. A new methodology for calculating precise and accurate values of E-s is presented. It is shown for the first time that E-s rather than E is constant during structural changes occurring within a supercooled liquid. It is also explained how to prepare and conduct microsecond range bulk-phase NVE simulations with essentially zero energy drift without the need for thermostating. The drift is analyzed with block averaging and new drift functions of the shadow energy. With such minimal drift, extremely small and accurate standard errors in the mean for quantities like E-s, E, and temperature, T, can be obtained. Values of the standard error for E-s of approximate to 10(-10) in molecule-based reduced units can be routinely achieved for simulations of 10(8) time steps. This corresponds to a simulation temperature drift of approximate to 10(-6) K/mu s, six orders of magnitude smaller than generally considered to be acceptable for protein simulations. We also show for the first time how these treatments can be extended with no loss of accuracy to polyatomic systems with both flexible degrees of freedom and arbitrary geometric constraints imposed via the SHAKE algorithm. As a bonus, estimates of simulation-average kinetic and total energies from high order velocity expressions can be obtained to a good approximation from 2nd order velocities and the average mean square force (for polyatomics, this refers to per site, including any constraint forces).