The previous notes sampled the ring-polymer distribution by molecular dynamics. Path-integral Monte Carlo targets the same configurational distribution with Metropolis moves instead of fictitious momenta and timesteps. This note writes the local bead update, explains where constrained PIMC enters, and checks an observable mean for a harmonic oscillator.

Target Distribution

After Trotter discretization, the one-dimensional canonical partition function can be written as a configurational path integral

\[ Z_P \propto \int dq_1\cdots dq_P\, \exp[-S_P(q)], \]

with dimensionless Euclidean action

\[ S_P(q)= \sum_{j=1}^{P} \left[ \frac{m}{2\beta_P\hbar^2}(q_{j+1}-q_j)^2 +\beta_P V(q_j) \right], \qquad \beta_P=\frac{\beta}{P}, \quad q_{P+1}=q_1 . \]

PIMD samples this distribution by dynamics on an extended Hamiltonian. PIMC samples it directly as a probability density proportional to \(e^{-S_P}\). There is no MD stability question, but there is a Monte Carlo mixing question.

Local Metropolis Move

The implementation follows the simple local-action structure used in the source notes and in the toy-model code. During one sweep, every bead is visited once in a random order. If bead \(i\) is proposed to move from \(q_i\) to \(q_i'\), only the two adjacent spring links and the bead potential change:

\[ S_i(q_i)= \frac{m}{2\beta_P\hbar^2} \left[(q_i-q_{i-1})^2+(q_{i+1}-q_i)^2\right] +\beta_P V(q_i). \]

The move is accepted with

\[ p_{\mathrm{acc}}=\min\left(1,\exp[-(S_i(q_i')-S_i(q_i))]\right). \]

For a position observable \(A(q)\), the estimator is evaluated on retained paths after burn-in:

\[ A_P(q)=\frac{1}{P}\sum_{j=1}^{P}A(q_j). \]

Constrained PIMC Context

The same sampler can be modified to keep the centroid \(q_c=P^{-1}\sum_j q_j\) fixed. A convenient proposal moves one bead and then shifts the whole path so that the centroid returns to a target value \(x_c\). This is the starting point for centroid effective-potential or CMD-style calculations:

\[ V_c(x_c;\beta)= -\frac{1}{\beta} \ln \left\langle \exp\left[-\beta_P\sum_{j=1}^{P}V(q_j)\right] \right\rangle_{\mathrm{free},\,q_c=x_c} +\mathrm{const}. \]

This Part IV keeps the benchmark unconstrained because the goal is to verify the basic equilibrium PIMC sampler first. The constrained formula is included to show how the method extends beyond ordinary thermal sampling.

Workflow

  1. Use the same harmonic oscillator scale as Part I: \(m=\hbar=1\), \(\beta=2\), \(\omega_0=4\), and \(P=32\).
  2. Start eight independent chains from \(q_j=0\).
  3. Run local single-bead Metropolis sweeps with adaptive proposal size during burn-in, targeting an acceptance rate near 0.45.
  4. After burn-in, retain every fifth path and evaluate the bead-average potential \(V_P(q)=P^{-1}\sum_j m\omega_0^2q_j^2/2\).
  5. Compare the Monte Carlo mean with both the exact finite-\(P\) harmonic value and the infinite-\(P\) quantum value.

Result

For \(P=32\), the run gives \(\langle V\rangle_{\mathrm{PIMC}}=0.99021\pm0.00684\). The exact finite-\(P\) target is \(0.99296\), while the infinite-\(P\) quantum value is \(1.00067\). Within the Monte Carlo uncertainty, the sampler is targeting the finite-bead ring-polymer distribution as intended.

PIMC harmonic oscillator Metropolis benchmark
Metropolis PIMC benchmark for the harmonic oscillator with \(m=\hbar=1\), \(\beta=2\), \(\omega_0=4\), and \(P=32\). Left: horizontal axis is the number of retained paths on a log scale; vertical axis is the running mean of the bead-average potential estimator \(\langle V\rangle\). The red curve is the PIMC running mean, the gray dashed line is the exact finite-\(P\) target \(0.99296\), the black line is the exact quantum value \(1.00067\), and the pale red band is the final PIMC mean plus or minus two batch-standard-error estimates. Right: histogram of all retained bead coordinates. The gray dashed curve is the finite-\(P\) Gaussian bead marginal, and the black curve is the infinite-\(P\) quantum Gaussian; for this \(P=32\) benchmark the two nearly overlap.

Interpretation

PIMC and PIMD share the same finite-\(P\) configurational target. PIMD adds fictitious momenta and then worries about thermostats, normal-mode stiffness, and free-step stability. PIMC removes the timestep issue and replaces it with proposal design, burn-in, acceptance-rate control, and correlation between retained paths. For a harmonic oscillator, the finite-\(P\) target is known analytically, so it is a clean test that the Metropolis action and observable estimator are implemented correctly.

Code used in this note

  • pimc_sho_metropolis.py implements the local-action Metropolis sampler, runs the SHO benchmark, writes the figure, and saves summary tables.
  • pimc-sho-summary.csv stores the PIMC mean, finite-\(P\) target, exact quantum value, \(\langle q^2\rangle\), acceptance rate, and final proposal size.
  • pimc-sho-running.csv stores the running \(\langle V\rangle\) curve used in the left panel.

References

  • Feynman and Hibbs, Quantum Mechanics and Path Integrals.
  • Ceperley, Rev. Mod. Phys. 67, 279 (1995).
  • Tuckerman, Statistical Mechanics: Theory and Molecular Simulation.

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.