Module 8 · Quantum Dynamics
Ring Polymer Molecular Dynamics
Ring Polymer Molecular Dynamics (RPMD) is a rigorous method for computing quantum mechanical rate constants from classical MD trajectories of an extended “ring polymer” system. It is the reference method for quantum tunneling corrections in this course, against which MAML-NequIP predictions are benchmarked.
8.1 The Ring Polymer Representation
Each quantum nucleus is replaced by a ring of \(n\) classical “beads” connected by harmonic springs, obtained from the Feynman path-integral representation of the quantum partition function:
\[ Z_n = \left(\frac{m n}{2\pi\beta\hbar^2}\right)^{nN/2} \int \prod_{\alpha=1}^n d\mathbf{R}^{(\alpha)}\, \exp\!\left(-\beta_n \sum_\alpha \left[\frac{mn}{2\beta^2\hbar^2}\left|\mathbf{R}^{(\alpha)}-\mathbf{R}^{(\alpha-1)}\right|^2 + V\!\left(\mathbf{R}^{(\alpha)}\right)\right]\right) \tag{8.1}\]
where \(\beta_n = \beta/n\), and the ring closure condition \(\mathbf{R}^{(n+1)} = \mathbf{R}^{(1)}\) is imposed. In the limit \(n \to \infty\), Eq. (8.1) becomes exact.
8.2 The RPMD Rate Constant
The RPMD rate is computed via the Bennett–Chandler two-step procedure:
\[ k_{\text{RPMD}} = k_{\text{TST}}^{\text{RPMD}} \cdot \kappa_{\text{RPMD}} \tag{8.2}\]
where \(k_{\text{TST}}^{\text{RPMD}}\) is the RPMD transition-state-theory rate (computed from the centroid PMF) and \(\kappa_{\text{RPMD}}\) is the recrossing transmission coefficient (from forward flux sampling).
8.3 Practical Setup
# i-PI RPMD input — key parameters
temperature: 300 K
n_beads: 32 # for H at 300K; use 48 at 250K
timestep: 0.25 fs
thermostat: PILE # Path Integral Langevin Equation
friction: 1/300 fs⁻¹
# Collective variable: proton transfer coordinate
xi: d(D-H) - d(H-A) # ξ ∈ [-0.75, 0.75] Å
# Umbrella sampling windows
n_windows: 30
delta_xi: 0.05 Å
k_spring: 50 kcal/mol/Ų
t_equil: 20 ps
t_prod: 100 ps
# Force calculator: MAML-NequIP via ASE interface
calculator: MAMLNequIPCalculator(adapted_model=phi_IFP)8.4 Number of Beads: Temperature Scaling
The required number of beads scales as \(n \sim \hbar\omega_{\max}/k_BT\), where \(\omega_{\max}\) is the highest-frequency mode of the transferring nucleus. For H transfer at 300 K: \(n \approx 32\)–64. For D, scale by \(\sqrt{m_H/m_D} \approx 0.71\): \(n \approx 24\)–48. The computational cost of RPMD scales linearly with \(n\) when the force calculator (MAML-NequIP) is the bottleneck.