Part V used normal modes to propagate an RPMD harmonic benchmark. This note steps back to the integrator itself. The i-PI normal-mode code diagonalizes the free ring-polymer spring matrix and propagates each mode analytically; the alternative bead-coordinate route integrates the same spring force by velocity Verlet, often with subcycling. The difference is small at low bead number, but it becomes the controlling stability issue when \(P\) is large.

Free Ring-Polymer Block

In a split PIMD or RPMD integrator, the stiff part is the free ring-polymer Hamiltonian

\[ H_P^0(q,p)= \sum_{j=1}^{P}\frac{p_j^2}{2m} +\frac{1}{2}m\omega_P^2 \sum_{j=1}^{P}(q_{j+1}-q_j)^2, \qquad \omega_P=\frac{P}{\beta\hbar}. \]

The spring matrix is circulant. Its Fourier eigenvectors diagonalize the internal motion and give

\[ \omega_k= 2\omega_P\sin\left(\frac{k\pi}{P}\right), \qquad k=0,\ldots,P-1. \]

In normal-mode coordinates, the free step is just independent harmonic motion. For \(\omega_k>0\),

\[ \begin{bmatrix} Q_k(t+\Delta t)\\ P_k(t+\Delta t) \end{bmatrix} = \begin{bmatrix} \cos\theta_k & \sin\theta_k/(m\omega_k)\\ -m\omega_k\sin\theta_k & \cos\theta_k \end{bmatrix} \begin{bmatrix} Q_k(t)\\ P_k(t) \end{bmatrix}, \qquad \theta_k=\omega_k\Delta t . \]

The centroid mode has \(\omega_0=0\) and drifts freely. This is the idea behind i-PI's `normalmodes.py`: transform to normal modes by FFT or by an explicit matrix, build the per-mode propagator, apply it, and transform back.

Bead-Coordinate Alternative

The same free spring block can be integrated directly in bead coordinates:

\[ F_j^{\mathrm{spring}}= -m\omega_P^2(2q_j-q_{j+1}-q_{j-1}). \]

A bead-coordinate velocity-Verlet substep is

\[ p\leftarrow p+\frac{\delta t}{2}F^{\mathrm{spring}}(q), \qquad q\leftarrow q+\delta t\,p/m, \qquad p\leftarrow p+\frac{\delta t}{2}F^{\mathrm{spring}}(q). \]

This is simple and avoids transforms. The price is that Verlet must resolve the highest spring frequency. Since \(\omega_{\max}\approx2\omega_P=2P/(\beta\hbar)\), the linear stability condition is roughly

\[ \delta t\,\omega_{\max}<2 \quad\Rightarrow\quad \delta t\lesssim\frac{\beta\hbar}{P}. \]

Therefore, at fixed outer timestep \(\Delta t\), the bead-coordinate route needs a subcycle count \(n_{\mathrm{mts}}\) that grows approximately linearly with \(P\). This is the role of the Cartesian `free_babstep` style path in the reference i-PI code.

Expected Limits

  • At \(P=1\), there are no internal modes, so normal-mode and bead-coordinate propagation are the same classical free drift.
  • At small \(P\), bead-coordinate Verlet can be cheaper because it avoids FFTs and a single step may still be stable.
  • At large \(P\), the largest spring frequency grows with \(P\). Fixed-step bead-coordinate Verlet either becomes inaccurate or unstable.
  • Subcycling restores stability, but its required \(n_{\mathrm{mts}}\) also grows with \(P\), so the apparent low-\(P\) efficiency advantage can disappear.
  • Normal-mode propagation removes the spring stiffness from the timestep restriction. It does not remove physical-force timestep limits, and Part III still matters for exact-free-step resonance and Cayley stabilization.

Workflow

  1. Use the free ring-polymer Hamiltonian only, with \(m=\hbar=1\), \(\beta=2\), outer timestep \(\Delta t=0.04\), and 600 outer steps.
  2. Initialize a deterministic high-mode-rich bead configuration and momentum for each bead count.
  3. Compare three propagators: FFT normal-mode exact propagation, bead-coordinate Verlet with \(n_{\mathrm{mts}}=1\), and bead-coordinate Verlet with automatic subcycling.
  4. Scan \(P=8,16,32,64,128,256,512\), record maximum relative free-ring-polymer energy error, and time the completed trajectories.
  5. At \(P=128\), scan \(n_{\mathrm{mts}}\) directly to see how bead-coordinate convergence and runtime change.

Result: Bead Number Scan

With \(\Delta t=0.04\), one-step bead-coordinate Verlet is stable through \(P=32\), but the first unstable case appears at \(P=64\). Normal-mode FFT propagation conserves the free-ring-polymer energy to roundoff; at \(P=512\), the maximum relative error is \(1.07\times10^{-13}\). The automatic bead-coordinate subcycling remains stable at high \(P\), but at \(P=512\) it uses \(n_{\mathrm{mts}}=26\), gives a maximum relative energy error of \(7.62\times10^{-2}\), and takes about \(0.407\) s for the fixed trajectory versus \(0.048\) s for the normal-mode route.

Normal-mode and bead-coordinate free ring-polymer propagation benchmark
Free ring-polymer propagation benchmark with \(m=\hbar=1\), \(\beta=2\), \(\Delta t=0.04\), and 600 outer steps. Left: horizontal axis is bead number \(P\); vertical axis is the maximum relative error in the conserved free-ring-polymer Hamiltonian. Black circles are FFT normal-mode exact propagation, red squares are bead-coordinate velocity Verlet with one outer step and no subcycling, and blue triangles are bead-coordinate Verlet with automatic \(n_{\mathrm{mts}}\) chosen from the highest spring frequency. The gray dotted line marks 100% relative energy error. Right: runtime for completed fixed-length trajectories. The red timing curve stops after \(P=32\) because the larger-\(P\) one-step bead trajectories blow up early and are not meaningful speed measurements.

Result: Subcycling Scan

At \(P=128\), bead-coordinate Verlet with \(n_{\mathrm{mts}}=1\) or 2 is unstable. The linear stability boundary predicts that \(n_{\mathrm{mts}}\ge3\) is needed for this timestep. Increasing \(n_{\mathrm{mts}}\) then reduces the energy error, but the cost grows almost linearly: the best tested case, \(n_{\mathrm{mts}}=24\), lowers the maximum relative energy error to \(6.58\times10^{-3}\), while remaining much slower than the exact normal-mode free step.

Subcycling convergence for bead-coordinate free ring-polymer propagation
Subcycling convergence at \(P=128\). Left: horizontal axis is bead-coordinate subcycle count \(n_{\mathrm{mts}}\); vertical axis is maximum relative free-ring-polymer energy error. The red curve is bead-coordinate velocity Verlet, the black dashed horizontal line is the FFT normal-mode energy error for the same trajectory, and the gray dotted vertical line is the approximate Verlet linear-stability boundary. Right: runtime for completed bead-coordinate trajectories as \(n_{\mathrm{mts}}\) increases; the black dashed line is the normal-mode runtime. Unstable bead-coordinate runs are omitted from the runtime panel.

Interpretation

The test mostly confirms the practical intuition, with one correction. At low bead count, bead-coordinate evolution can indeed be cheaper, and for \(P=8\) or \(16\) its errors may be acceptable for a loose free-step tolerance. But the advantage is conditional: as \(P\) grows, the spring frequency grows with it, so fixed-step bead-coordinate Verlet first becomes inaccurate and then unstable. Subcycling can recover robustness, but it turns the cost into a stiffness-management problem. Normal modes are the cleaner limit: the free spring block is solved analytically, and the remaining timestep questions come from the physical potential, thermostat choices, and the resonance issue discussed in Part III.

Code used in this note

References

  • Ceriotti, More, and Manolopoulos, Comput. Phys. Commun. 185, 1019 (2014).
  • Tuckerman, Berne, Martyna, and Klein, J. Chem. Phys. 99, 2796 (1993).
  • Korol, Bou-Rabee, and Miller, J. Chem. Phys. 151, 124103 (2019).

PIMD series map

  1. Part I: Ring-Polymer Basics - path-integral isomorphism, normal modes, estimators, and the SHO convergence test.
  2. Part II: NVT Thermostats and GLE - PILE, GLE, PIGLET-style covariance shaping, and sampling diagnostics.
  3. Part III: Cayley Stability - why exact free ring-polymer propagation can be unstable and how the Cayley step fixes it.
  4. Part IV: Path-Integral Monte Carlo - local Metropolis bead moves, constrained-PIMC context, and an SHO observable mean.
  5. Part V: From PIMD to RPMD - PIMC-initialized RPMD trajectories and an exact SHO \(q\)-\(q\) Kubo benchmark.
  6. Part VI: Normal-Mode Free Propagation - normal-mode and bead-coordinate free spring propagation, stability limits, and subcycling costs.