Path-integral molecular dynamics turns an equilibrium quantum Boltzmann problem into a classical sampling problem over a ring polymer. This first note keeps the scope deliberately static: it derives the ring-polymer Hamiltonian, writes the sampling equations of motion, and tests the bead-number convergence of a harmonic-oscillator potential estimator.

Quantum Boltzmann trace

For a one-dimensional Hamiltonian

\[ \hat H=\frac{\hat p^2}{2m}+V(\hat q), \]

the canonical partition function is

\[ Z=\operatorname{Tr}\left[e^{-\beta \hat H}\right] =\int dq\,\langle q|e^{-\beta\hat H}|q\rangle . \]

For a short imaginary-time interval \(\Delta\tau\), a symmetric kinetic–potential splitting gives a free-particle Gaussian kernel multiplied by two half-step potential weights:

\[ \begin{aligned} K_{\Delta\tau}(q,q')&\simeq \sqrt{\frac{m}{2\pi\hbar\Delta\tau}}\, e^{-g_{\Delta\tau}(q,q')},\\ g_{\Delta\tau}(q,q')&= \frac{m(q-q')^2}{2\hbar\Delta\tau}\\ &\quad+\frac{\Delta\tau[V(q)+V(q')]}{2\hbar}. \end{aligned} \]

Split the full Boltzmann operator into \(P\) factors with \(\Delta\tau=\beta\hbar/P\), insert a coordinate resolution between each pair, and take the trace. The trace identifies the final and initial coordinates, closing the path. Multiplication combines adjacent half-step potential factors and produces the spring-linked weight below; the omitted normalization is \([m/(2\pi\hbar\Delta\tau)]^{P/2}\).

Writing \(\beta_P=\beta/P\) and \(q_{P+1}=q_1\), the result is

\[ Z_P \propto \int dq_1\cdots dq_P\, \exp\left[ -\beta_P\sum_{j=1}^{P} \left\{ \frac{1}{2}m\omega_P^2(q_{j+1}-q_j)^2+V(q_j) \right\} \right], \qquad \omega_P=\frac{1}{\beta_P\hbar}. \]

The quantum particle is therefore represented by \(P\) classical beads connected by harmonic springs. Increasing \(P\) resolves shorter imaginary-time structure, but also introduces high-frequency internal modes.

Ring-polymer Hamiltonian

To sample the configurational distribution with molecular dynamics, fictitious bead momenta are introduced:

\[ H_P(q,p)= \sum_{j=1}^{P} \left[ \frac{p_j^2}{2m'} +\frac{1}{2}m\omega_P^2(q_{j+1}-q_j)^2 +V(q_j) \right]. \]

The mass \(m'\) and the MD time are sampling devices; they do not define the real-time quantum motion. The corresponding sampling equations are

\[ \dot q_j=\frac{p_j}{m'}, \qquad \dot p_j= -m\omega_P^2(2q_j-q_{j+1}-q_{j-1}) -\frac{\partial V(q_j)}{\partial q_j}. \]

In practice these equations are usually combined with thermostats, because the free ring-polymer modes are stiff and can be poorly ergodic. Part II focuses on that NVT layer.

Normal modes

The spring matrix is diagonal in ring-polymer normal modes. For the free ring polymer, the mode frequencies are

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

The \(k=0\) mode is the centroid. The remaining modes describe internal imaginary-time fluctuations. Splitting algorithms exploit this diagonal form by treating the free ring-polymer substep analytically or with a stable approximation.

Normal-mode EOM workflow

A practical PIMD trajectory usually does not integrate the stiff spring force directly in bead coordinates. One standard second-order step splits the physical potential from the free ring-polymer Hamiltonian:

  1. Apply a half kick from the physical or external potential in bead coordinates, \(p_j\leftarrow p_j-\frac{\Delta t}{2}\partial V(q_j)/\partial q_j\), using the current bead positions.
  2. Transform bead positions and momenta to normal modes, \(Q=Cq\) and \(P=Cp\), using an orthogonal ring-polymer transform.
  3. Propagate each free normal mode analytically. For \(\omega_k\ne0\),
\[ \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 . \]

For the centroid mode, where \(\omega_k=0\), this reduces to the free drift \(Q_0\leftarrow Q_0+\Delta t\,P_0/m'\).

  1. Transform the updated normal modes back to bead coordinates, \(q=C^TQ\) and \(p=C^TP\).
  2. Apply the second half kick from the physical or external potential, now evaluated at the updated bead positions: \(p_j\leftarrow p_j-\frac{\Delta t}{2}\partial V(q_j)/\partial q_j\).
  3. For NVT sampling, apply the thermostat step in the chosen splitting order, often in normal-mode coordinates; then accumulate estimators only after equilibration.

Thus the deterministic core is a symmetric external-force half step, full free-ring-polymer normal-mode step, and external-force half step. The benchmark below uses the same normal-mode structure but evaluates the harmonic equilibrium covariance directly, so it has no timestep error.

Static estimators

For a position-only observable \(O(q)\), the standard bead estimator is

\[ O_P(q)=\frac{1}{P}\sum_{j=1}^{P}O(q_j), \qquad \langle \hat O\rangle = \lim_{P\to\infty} \left\langle O_P(q)\right\rangle_{H_P}. \]

For a harmonic oscillator \(V(q)=m\omega_0^2q^2/2\), the exact quantum potential energy is

\[ \langle V\rangle_{\mathrm{exact}} = \frac{\hbar\omega_0}{4} \coth\left(\frac{\beta\hbar\omega_0}{2}\right). \]

The finite-\(P\) ring-polymer result can also be evaluated directly in normal modes:

\[ \langle V_P\rangle = \frac{1}{2\beta} \sum_{k=0}^{P-1} \frac{\omega_0^2}{\omega_k^2+\omega_0^2}. \]

Workflow

  1. Set \(m=\hbar=1\), \(\beta=2\), and \(\omega_0=4\).
  2. Compute the ring-polymer mode frequencies and the exact quantum value of \(\langle V\rangle\).
  3. Evaluate the finite-\(P\) normal-mode formula for \(P=1,\ldots,128\).
  4. Draw independent normal-mode Gaussian samples from the finite-\(P\) equilibrium distribution at powers of two.
  5. Scan all integer \(P\) values to estimate how many beads are needed for preset relative-error tolerances.
  6. Compare the estimator against the exact quantum value and the classical \(P=1\) limit.

Result

For this low-temperature oscillator, \(\beta\hbar\omega_0=8\), the exact finite-\(P\) formula gives a useful bead-count rule of thumb: the first integer \(P\) below 1% relative error is \(P=29\), the first below 0.2% is \(P=64\), and the first below 0.1% is \(P=90\). If only powers of two are used, this means \(P=32\), \(64\), and \(128\), respectively.

PIMD harmonic oscillator potential estimator convergence
Harmonic-oscillator PIMD benchmark. Left: horizontal axis is bead number \(P\) on a base-2 log scale; vertical axis is \(\langle V\rangle\). The black line is the exact quantum value \(1.000671\), the gray dashed line is the classical \(P=1\) value \(0.25\), the red curve is the finite-\(P\) normal-mode formula, and blue squares are direct normal-mode samples with two-standard-error bars. Right: relative error of the finite-\(P\) estimator. The dotted gray guides mark analytic convergence thresholds: 1% error at \(P=29\) and 0.2% error at \(P=64\).

Code used in this note

References

  • Feynman and Hibbs, Quantum Mechanics and Path Integrals.
  • Tuckerman, Berne, Martyna, and Klein, J. Chem. Phys. 99, 2796 (1993).
  • Ceriotti, Parrinello, Markland, and Manolopoulos, J. Chem. Phys. 133, 124104 (2010).

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.