VESPER: Online Learning of Free Energies with a Multi-Tempered Variational Objective
Authors: Jay L. Kaplan, Gustavo R. Pérez-Lemus, Mengjian Hua, Juan J. de Pablo, Eric Vanden-Eijnden
Organizations: Courant Institute School of Mathematics, Computing, and Data Science, New York University, New York, NY 10012, USA · Department of Chemical and Biomolecular Engineering, Tandon School of Engineering, New York University, Brooklyn, NY 11201, USA · Department of Physics, New York University, New York, NY 10012, USA · Machine Learning Lab, Capital Fund Management
Computing molecular free-energy landscapes over collective variables (CVs) is central to understanding molecular mechanisms, but poses a difficult learning problem without training data available a priori. Samples of the marginal CV distribution must be generated by molecular dynamics (MD), whose trajectories can remain trapped in metastable states and rarely cross the free-energy barriers governing conformational changes or chemical reactions. We introduce VESPER, which couples the generation of these samples to neural-network-driven free-energy learning: MD supplies configurations for learning, while the evolving model biases the MD to accelerate exploration and generate further training data. VESPER minimizes a single variational objective that is convex in the free energy and can combine several tempering levels, each flattening the free-energy landscape to promote barrier crossing. Every level has the same physical free energy as its minimizer, allowing data from all levels to train one shared model without reweighting. Learning uses the difference in mean parameter gradients between biased MD configurations and auxiliary walkers that explore the flattened model in CV space. We prove conditional convergence of the coupled overdamped population dynamics without replica exchange under log-Sobolev and Polyak-Lojasiewicz assumptions. VESPER matches or improves on established enhanced-sampling methods on alanine dipeptide and yields self-consistent free-energy reconstructions from three to six CVs on alanine tetrapeptide. For methane dehydrogenation on Ni(511) with a machine-learned interatomic potential, it recovers two-CV Boltzmann marginals and a reactant-to-saddle free-energy difference along a three-CV path, both qualitatively consistent with the lower-dimensional reference. It also samples RMSD-defined folded and unfolded regions of chignolin while biasing seventeen CVs.
Figures & tables
Figure 1: VESPER couples sampling and learning. Each tempering level k carries two walker populations, and all levels share one model Gθ . Left: physical walkers evolve in configuration space under U(x)−(1−αk)Gθ(s(x)) and are mapped to CV space by s(x) . Replica exchange swaps configurations between adjacent levels. Right: auxiliary walkers evolve directly in CV space under αkGθ(s) . The two populations supply the averages in equation 10 , whose ck -weighted difference updates θ . The updated model then sets the forces on both populations. At convergence, both populations at level k sample the same flattened landscape αkG⋆ . During learning, the instantaneous target for the physical CVs is QGθ,k∝pe(β−β^k)Gθ .
Figure 2: Analytic landscapes. (A) Full-evaluation-domain RMSE in kBT and (B) dimensionless DKL(ρ∥ρθ) against reduced time on Müller–Brown. Curves show means for single levels and a five-level ladder with (L-RE) and without (L-NRE) exchange, with 1,000 walkers per population. (C) Reference and learned surfaces at t=80 on a common 0 – 20kBT scale. Red and blue dots are physical and auxiliary walkers. (D–F) Four-well landscape in d=2 – 100 coordinates, with fixed hidden width and primary walker budget: (D) final SWA full-coordinate NRMSE, (E) DKL(ρ∥ρθ) , and (F) reduced time to sustained EMA NRMSE below 0.25 . Points and error bars show means and sample SDs. Open circles in (D,E) mark values above the plotted range or non-finite divergences. Open markers in (F) indicate at least one run without a sustained crossing.
Figure 3: Alanine dipeptide. (A) RMSE against the reference surface [ Zubieta Rico et al., 2025 ] over 0–10 ns per replica, averaged over three independent runs. VESPER curves are solid and baselines dashed. Gray shading marks errors above the chemical accuracy of 1 kcal mol -1 (4.184 kJ mol -1 ). L-RE, L-NRE: four-level ladders with and without exchange; OPES: on-the-fly probability enhanced sampling [ Invernizzi & Parrinello, 2020 ] ; MetaD: metadynamics [ Laio & Parrinello, 2002 ] ; VES: neural variationally enhanced sampling [ Valsson & Parrinello, 2014 , Bonati et al., 2019 ] ; TAMD: temperature-accelerated MD [ Maragliano & Vanden-Eijnden, 2006 , Bajpai et al., 2023 ] ; CFF: combined force-frequency sampling [ Sevgen et al., 2020 ] ; S-CFF: its SIREN variant [ Zubieta Rico et al., 2025 ] . (B) Reference (top left), and VESPER surfaces for α=1 (top right), a single level at α=0.01 (bottom left) and the exchange ladder (bottom right), on a common scale. Red and blue dots are physical and auxiliary walkers. (C) ADP with ϕ,ψ labeled [ Humphrey et al., 1996 ] . (D) ϕ of the physical walkers over 1–10 ns, colored by the learned free energy.
Figure 4: Alanine tetrapeptide: consistency across CV dimension. Two-dimensional Boltzmann marginals of (A) the three-CV and (B) the six-CV model, each from the mean surface of the four-level ladder with exchange. Columns: (ϕ1,ϕ2) , (ϕ1,ϕ3) , (ϕ2,ϕ3) . (C) ATP with ϕ bonds in blue and ψ bonds in red.
Figure 5: Methane on Ni(511): three-CV free energy. (A–C) Two-dimensional Boltzmann marginals in the planes (C–Ni, C–H coordination), ( S1 , C–H coordination) and ( S1 , C–Ni coordination). White line: bounded path on the mean surface obtained by a string-method calculation [ E et al., 2007 ] followed by constrained gradient-flow refinement, from reactant (1) through the saddle (3) to product (5). (D) Trajectory frames corresponding to the labeled points.
Figure 6: Chignolin with seventeen dihedral CVs (four independent runs, twelve levels from α=1 to 1/18 ). (A) Mean number of levels holding folded (C α RMSD below 0.20 nm), intermediate (labeled Partial) or unfolded (above 0.40 nm) configurations. Shading shows the minimum-to-maximum range across runs. (B) Mean folded fraction at each tempering level. In A and B, occupancies are averaged in 0.2 ns bins and across the four runs. Both panels display the first 25 ns per replica. (C) Free-energy projection at the end of the run from unbiased-level auxiliary walkers reconstructed on a fixed backbone, in C α radius of gyration and RMSD (dots: MD frames, red, and auxiliary walkers, blue, from all levels). (D) Dihedral marginal in (ψ3,ψ4) from the same unbiased-level auxiliary walkers. Contours every 5 kJ mol -1 .
Appendix figures & tables6 assets
Supplementary material from the paper’s appendix.
Appendix
System
α
K
Physical
Auxiliary
ADP
0.1,0.05,0.02,0.01
4
1
104
ATP, D=3,6
0.1,0.05,0.02,0.01
4
1
2×104
Methane/Ni(511)
0.05
1
1
6.5×104
Chignolin
18−k/11 , k=0,…,11
12
1
105
Appendix
Table 1: Levels per molecular system: tempering exponents α , number of levels K , and walkers per level.
System
Method
K
ns/replica
aggregate ns
Runs
ADP
VESPER, α=0.01 or 1
1
10
10
3
ADP
VESPER ladders (L-RE, L-NRE)
4
10
40
3
ADP
Each baseline
1
10
10
3
ATP, D=3,6
VESPER, α=0.01
1
10
10
3
ATP, D=3,6
VESPER ladders (L-RE, L-NRE)
4
10
40
3
Methane/Ni(511)
VESPER, α=0.05
1
7.5
7.5
3
Appendix
Table 2: Replicas K , simulated time per replica at the readout, and aggregate MD per run ( K times the time per replica).
System
Quantity
Value
ADP
RMSE, VESPER α=1 (kJ mol -1 )
13.27±0.26
ADP
RMSE, VESPER α=0.01
2.19±0.02
ADP
RMSE, L-RE / L-NRE
2.02±0.24 / 2.05±0.09
ADP
RMSE, MetaD/VES/S-CFF
2.61±0.22 / 4.72±0.25 / 2.04±0.51
ATP
r2 , 3-CV vs 6-CV marginals, three planes
0.951 / 0.895 / 0.890
ATP
RMSE of the same planes (kJ mol -1 )
5.5 / 8.1 / 7.5
Appendix
Table 3: Values reported in the main text for the molecular systems. ADP: 10 ns per replica, SWA8 readout, three independent runs. ATP: 10 ns, SWA10, three independent runs, mean surfaces for the plane comparisons. Methane: 7.5 ns, EMA readout, three independent runs. Its 1.12±0.25 eV entry is the route value on the mean surface with the sample SD evaluated at fixed route points. Chignolin: event counts use the first 25 ns per replica and the unbiased-level occupancy uses the full 26.2 ns trajectory, four independent runs. RMSEs are mean ± sample SD over runs.
α
RMSE
Walkers
RMSE
0.05
0.109±0.018
1
1.948±0.749
0.1
0.098±0.026
10
0.377±0.160
0.25
0.104±0.024
100
0.169±0.082
0.5
0.093±0.029
1000
0.104±0.024
1
0.602±0.494
10000
0.070±0.017
Appendix
Table 4: Müller–Brown errors at t=320 : Boltzmann-weighted FES RMSE in kBT of the parameter average, mean ± sample SD over four independent runs, against the tempering exponent at 1000 walkers per population and against the walker count at α=0.25 . A flat surface scores 1.281 .
NRMSE
DKL
d
L-RE
Best single ( α )
L-NRE
L-RE
Best single ( α )
L-NRE
2
0.36±0.20
0.28±0.10(0.05)
0.38±0.11
0.07±0.07
0.05±0.03(0.05)
0.07±0.03
5
0.42±0.10
0.26±0.09(0.5)
0.41±0.31
0.18±0.07
0.08±0.05(0.5)
0.22±0.27
10
0.23±0.07
0.16±0.04(0.25)
0.22±0.06
0.13±0.07
0.06±0.02(0.25)
0.12±0.06
25
0.22±0.06
0.20±0.03(0.5)
0.16±0.01
0.27±0.13
0.24±0.08(0.5)
0.18±0.01
50
0.19±0.01
0.24±0.04(0.5)
0.17±0.01
0.41±0.04
0.60±0.09(0.5)
0.36±0.04
Appendix
Table 5: Four-well dimensional scaling: mean ± sample SD of the full-coordinate NRMSE and DKL(ρ∥ρθ) at the final SWA readout. Each best single level is selected by the lowest mean of the corresponding metric among the tested levels at that dimension. All primary arms use the same total walker budget.
NRMSE
DKL
d
L-RE
L-NRE
L-RE
L-NRE
2
0.14±0.03
0.18±0.06
0.010±0.004
0.02±0.01
5
0.26±0.02
0.24±0.04
0.08±0.01
0.07±0.02
10
0.27±0.03
0.26±0.02
0.16±0.03
0.17±0.03
25
0.25±0.05
0.20±0.05
0.40±0.20
0.24±0.13
50
0.20±0.02
0.179±0.001
0.47±0.09
0.39±0.01
Appendix
Table 6: Larger-population controls for the four-well landscape: mean ± sample SD with 1,000 walkers per level and population, four times the primary ladder budget.
Molecular systems have many degrees of freedom, but their metastable behavior can often be described by a few collective variables. Identifying these variables and estimating free energies along them from limited simulation data remains a challenging, important problem. Separate short trajectories may sample different metastable states without capturing transitions or establishing their relative equilibrium populations. For unbiased trajectories generated with the same Hamiltonian at a single temperature, alternate methods based on histogram reweighting cannot correct this imbalance. Here we present SPIBER, which combines the State Predictive Information Bottleneck (SPIB) with Generative Flow Networks (GFlowNets). SPIB uses deep learning to approximate slow degrees of freedom through a past-future information bottleneck, retaining information needed to predict future metastable states. We show that this compression limits conditional entropy variations in populated regions, allowing conditional mean potential energies, which are much easier to calculate, to be used to approximate free energy differences. Given sufficient local sampling to estimate these energies, they define the target distribution for GFlowNets, energy-based generative samplers that sample according to estimated thermodynamic stability rather than observed populations. For a particle in a radial double-well potential, for alanine dipeptide, and for the nine-residue peptide AIB9, SPIBER recovers free energy differences between sampled metastable states to within one thermal energy unit of reference values. The method combines collective-variable learning and free energy estimation in up to four latent dimensions, without requiring converged state populations or additional molecular dynamics simulations.
Venkata Sai Sreyas Adury, Pratyush Tiwary
Chemical Physics Program and Institute for Physical Science and Technology, University of Maryland, College Park 20742, USA · Biophysics Program and Institute for Physical Science and Technology, University of Maryland, College Park 20742, USA · Department of Chemistry and Biochemistry and Institute for Physical Science and Technology, University of Maryland, College Park 20742, USA +1
Sampling the free energy surface, namely, the distribution of collective variables (CVs), is a crucial problem in statistical physics, as it underpins a better understanding of chemical reactions and conformational transitions. Traditional methods for free energy surface sampling involve simulation in high-dimensional configuration space and projecting the resulting configurations onto the CV space. To reduce the computational costs of such sampling, we propose FES-FM, a reduced flow matching (FM) method for free energy sampling (FES). We train a dynamical transport map in the CV space, thereby enabling direct sampling of the free energy surface. For many-particle systems, we construct a prior distribution based on the Hessian at a local minimum of the potential, which ensures both rotation-translation invariance and physically meaningful configurations. We evaluate the proposed method across a variety of potential functions and collective variables. Comparative experiments demonstrate that our approach drastically reduces computational costs while delivering superior accuracy per unit sampling time.
Zichen Liu, Tiejun Li
Center for Data Science, Peking University · LMAM and School of Mathematical Sciences, Peking University · Center for Machine Learning Research, Peking University
Rare events are central to the evolution of complex many-body systems, characterized as key transitional configurations on the free energy surface (FES). Conventional methods require adequate sampling of rare event transitions to obtain the FES, which is computationally very demanding. Here we introduce the variational free energy surface (VaFES), a dataset-free framework that directly models FESs using tractable-density generative models. Rare events can then be immediately identified from the FES with their configurations generated directly via one-shot sampling of generative models. By extending a coarse-grained collective variable (CV) into its reversible equivalent, VaFES constructs a latent space of intermediate representation in which the CVs explicitly occupy a subset of dimensions. This latent-space construction preserves the physical interpretability and transparent controllability of the CVs by design, while accommodating arbitrary CV formulations. The reversibility makes the system energy exactly accessible, enabling variational optimization of the FES without pre-generated simulation data. A single optimization yields a continuous, differentiable FES together with one-shot generation of rare-event configurations. Our method can reproduce the exact analytical solution for the bistable dimer potential and identify a chignolin native folded state in close alignment with the experimental NMR structure. Our approach thus establishes a scalable, systematic framework for advancing the study of complex statistical systems.
Shuo-Hui Li, Chen Chen, Yao-Wen Zhang +1
Department of Physics, The Hong Kong University of Science and Technology, Hong Kong, China · Department of Chemistry, The Hong Kong University of Science and Technology, Hong Kong, China · IAS Center for AI for Scientific Discoveries, Hong Kong University of Science and Technology, Hong Kong, China