  Path-integral-based molecular dynamics (MD) simulations are widely used for
the calculation of numerically exact quantum Boltzmann properties and
approximate dynamical quantities. A nearly universal feature of MD numerical
integration schemes for equations of motion based on imaginary-time path
integrals is the use of harmonic normal modes for the exact evolution of the
free ring-polymer positions and momenta. In this work, we demonstrate that this
standard practice creates numerical artifacts. In the context of conservative
(i.e., microcanonical) equations of motion, it leads to numerical instability.
In the context of thermostatted (i.e., canonical) equations of motion, it leads
to non-ergodicity of the sampling. These pathologies are generally proven to
arise at integration timesteps that depend only on the system temperature and
the number of ring-polymer beads, and they are numerically demonstrated for the
cases of conventional ring-polymer molecular dynamics (RPMD) and thermostatted
RPMD (TRPMD). Furthermore, it is demonstrated that these numerical artifacts
are removed via replacement of the exact free ring-polymer evolution with a
second-order approximation based on the Cayley transform. The Cayley
modification introduced here can immediately be employed with almost every
existing integration scheme for path-integral-based molecular dynamics -
including path-integral MD (PIMD), RPMD, TRPMD, and centroid MD - providing
strong symplectic stability and ergodicity to the numerical integration, at no
penalty in terms of computational cost, algorithmic complexity, or accuracy of
the overall MD timestep. Furthermore, it is shown that the improved numerical
stability of the Cayley modification allows for the use of larger MD timesteps.
We suspect that the Cayley modification will therefore find useful application
in many future path-integral-based MD simulations.
