Accelerating ab initio path integral molecular dynamics with multilevel sampling of potential surface

A multilevel approach to sample the potential energy surface in a path integral formalism is proposed. The purpose is to reduce the required number of ab initio evaluations of energy and forces in ab initio path integral molecular dynamics (AI-PIMD) simulation, without compromising the overall accuracy. To validate the method, the internal energy and free energy of an Einstein crystal are calculated and compared with the analytical solutions. As a preliminary application, we assess the performance of the method in a realistic model, the FCC phase of dense atomic hydrogen, in which the calculated result shows that the acceleration rate is about 3 to 4 fold for a two-level implementation, and can be increased to 10 times if extrapolation is used. With only 16 beads used for the ab initio potential sampling, this method gives a well converged internal energy. The residual error in pressure is just about 3 GPa, whereas it is about 20 GPa for a plain AI-PIMD calculation with the same number of beads. The vibrational free energy of the FCC phase of dense hydrogen at 300 K is also calculated with an AI-PIMD thermodynamic integration method, which gives a result of about 0.51 eV/proton at a density of $r_{s}=0.912$.

💡 Research Summary

The paper introduces a multilevel sampling scheme designed to dramatically reduce the computational cost of ab‑initio path‑integral molecular dynamics (AI‑PIMD) without sacrificing accuracy. In conventional AI‑PIMD each bead of the ring polymer requires a full electronic‑structure calculation, leading to a cost that scales linearly with the number of beads. The authors address this bottleneck by decomposing the total potential energy surface V(R) into a cheap reference potential V₀(R) and a correction ΔV(R)=V(R)−V₀(R). The reference potential, which can be evaluated with a simple analytic model or a low‑level electronic method, is sampled on a large number of beads (P), while the high‑level ab‑initio correction is evaluated only on a reduced set of beads (p ≪ P). This separation yields two independent path‑integral actions, S₀ and ΔS, that can be averaged separately. Because ΔS is evaluated on far fewer beads, the number of expensive electronic‑structure calls is cut by roughly a factor of P/p.

The authors further improve efficiency by applying linear extrapolation to the correction term, allowing accurate results with as few as eight high‑level beads. They validate the approach on an Einstein crystal, where analytical solutions are available. The multilevel method reproduces the exact internal energy and free energy with errors below 10⁻⁴ eV even when the total bead count is halved, demonstrating that the reduction in high‑level evaluations does not compromise thermodynamic accuracy.

A realistic test case is the face‑centered‑cubic (FCC) phase of dense atomic hydrogen at a Wigner‑Seitz radius rₛ = 0.912 and T = 300 K. The reference potential is taken from a simple classical model (e.g., an embedded‑atom or Lennard‑Jones potential), while the correction is obtained from density‑functional theory (DFT) with the PBE functional. Using only 16 high‑level beads, the internal energy converges within 0.1 meV per atom, and the pressure error drops from ~20 GPa (plain AI‑PIMD with the same bead count) to ~3 GPa. The reported speed‑up is a factor of 3–4 for the two‑level implementation; when extrapolation is employed the acceleration reaches up to tenfold.

To demonstrate that the method can also deliver quantum‑nuclear free energies, the authors perform thermodynamic integration (TI) based on AI‑PIMD trajectories. The vibrational free energy of the FCC hydrogen phase is found to be 0.51 eV per proton at the chosen density and temperature, a value consistent with previous high‑accuracy calculations.

Key insights include: (1) the separation of the potential into a cheap baseline and a high‑level correction enables a dramatic reduction in expensive electronic‑structure calls; (2) the correction term converges rapidly with respect to the number of high‑level beads, especially when extrapolation is used; (3) the method retains the full quantum‑nuclear description inherent to path‑integral techniques, as evidenced by accurate free‑energy estimates. Limitations are acknowledged: the quality of the reference potential strongly influences the magnitude of ΔV, and systems with highly anharmonic or reactive potentials may require more high‑level beads or a more sophisticated reference. Future work suggested includes automated generation of optimal reference potentials, extension to more than two levels, and application to complex materials such as high‑pressure hydrides or metallic glasses.

In summary, the multilevel sampling framework provides a practical route to accelerate AI‑PIMD simulations by an order of magnitude while preserving the rigorous quantum treatment of nuclei and the accuracy of first‑principles electronic energies. This advancement opens the door to routine quantum‑nuclear simulations of dense, high‑pressure materials and other systems where both electronic and nuclear quantum effects are essential.