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 ≈10^-10 in molecule-based reduced units can be routinely achieved for simulations of 10⁸ time steps. This corresponds to a simulation temperature drift of ≈10^-6 K/μ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).