Deriving Reliable Nucleation Rates from Metadynamics Simulations: Application to Yukawa Fluids

In order to solidify the usefulness of metadynamics in studying nucleation of crystals from supercooled liquids, we provide a specific procedure to calculate nucleation free energy barriers. After a pedagogical review of the important elements of classical nucleation theory and how metadynamics is used to find nucleation free energy barriers, we explain the benefits of local collective variables over more common global collective variables. We show how a metadynamics free energy barrier must be carefully postprocessed so that classical nucleation theory can be applied to calculate nucleation rates. We apply our procedure to a Yukawa plasma and show that a particular physically-motivated fit to metadynamics data reproduces low-temperature reference data, justifying the usefulness of metadynamics to predict nucleation rates and the nucleation critical temperature.

💡 Research Summary

The paper establishes a rigorous protocol for extracting reliable crystal nucleation rates from metadynamics simulations and demonstrates its effectiveness on a Yukawa one‑component plasma (YOCP). After a concise review of classical nucleation theory (CNT) and its limitations, the authors argue that the key to quantitative predictions lies in (i) choosing appropriate collective variables, (ii) properly reweighting the biased metadynamics data, and (iii) fitting “effective” CNT parameters (Δμ and γ) to the free‑energy surface obtained from metadynamics.

First, the authors introduce local collective variables that quantify the crystalline environment of each particle rather than relying on a global order parameter. This choice captures the often non‑spherical, diffuse, and defect‑rich nature of early nuclei, which global variables tend to miss. The metadynamics bias potential Vb is built as a sum of Gaussians deposited during the simulation; because Vb evolves with time, a simple subtraction of Vb from the biased free energy would be inaccurate. The authors therefore employ a balanced exponential reweighting scheme that accounts for the time‑dependence of Vb, yielding an unbiased free‑energy profile F(N) as a function of the nucleus size N (the number of solid‑like particles).

With the unbiased F(N) in hand, they fit it to the CNT functional form ΔF(N)=Δμ N+γ N2/3. The fitted Δμ and γ are termed “effective” because they incorporate deviations from the ideal CNT assumptions (e.g., curvature corrections, interfacial roughness, internal stresses). From the fitted curve they extract the critical nucleus size N* and the barrier height ΔF(N*). The kinetic prefactor D⁺, representing the rate of particle attachment/detachment at the cluster surface, is obtained from separate equilibrium simulations. Inserting these quantities into the standard CNT rate expression J = D⁺ s |Δμ|/(6πkBT N*) exp