Abstract
Convergence with respect to imaginary-time discretization (i.e., the number of ring-polymer beads) is an essential part of any path-integral-based molecular dynamics (MD) calculation. However, an unfortunate property of existing non-preconditioned numerical integration schemes for path-integral molecular dynamics-including essentially all existing ring-polymer molecular dynamics (RPMD) and thermostatted RPMD (T-RPMD) methods-is that for a given MD time step, the overlap between the exact ring-polymer Boltzmann-Gibbs distribution and that sampled using MD becomes zero in the infinite-bead limit. This has clear implications for hybrid Metropolis Monte Carlo/MD sampling schemes, and it also causes the divergence with bead number of the primitive path-integral kinetic-energy expectation value when using standard RPMD or T-RPMD. We show that these and other problems can be avoided through the introduction of "dimension-free" numerical integration schemes for which the sampled ring-polymer position distribution has non-zero overlap with the exact distribution in the infinite-bead limit for the case of a harmonic potential. Most notably, we introduce the BCOCB integration scheme, which achieves dimension freedom via a particular symmetric splitting of the integration time step and a novel implementation of the Cayley modification [R. Korol et al., J. Chem. Phys. 151, 124103 (2019)] for the free ring-polymer half-steps. More generally, we show that dimension freedom can be achieved via mollification of the forces from the external physical potential. The dimension-free path-integral numerical integration schemes introduced here yield finite error bounds for a given MD time step, even as the number of beads is taken to infinity; these conclusions are proven for the case of a harmonic potential and borne out numerically for anharmonic systems that include liquid water. The numerical results for BCOCB are particularly striking, allowing for nearly three-fold increases in the stable time step for liquid water with respect to the Bussi-Parrinello (OBABO) and Leimkuhler (BAOAB) integrators, while introducing negligible errors in the calculated statistical properties and absorption spectrum. Importantly, the dimension-free, non-preconditioned integration schemes introduced here preserve ergodicity and global second-order accuracy, and they remain simple, black-box methods that avoid additional computational costs, tunable parameters, or system-specific implementations.