Collective-Variable Stein Variational Gradient Descent for Free Energy Surface Exploration

demonstrate the potential of CV-SVGD as a flexible, parallelizable framework for enhanced sampling in multidimensional CV spaces. We additionally provide an open-source Python implementation to facilitate its adoption.

Introduction

Free-energy calculations are central to the use of molecular simulations for studying conformational changes, molecular recognition, chemical reactions, and other processes governed by rare events.1,2 In principle, equilibrium properties can be obtained by sampling the Boltzmann distribution using molecular dynamics (MD) or Monte Carlo simulations. In practice, however, transitions between metastable states may occur on timescales far beyond those accessible to straightforward simulations, motivating the development of a broad range of enhanced-sampling methods.3–48

A common strategy underlying many enhanced-sampling approaches is to describe the molecular system using a reduced set of collective variables (CVs) that capture the slow or functionally relevant degrees of freedom. When an appropriate low-dimensional representation is available, methods such as umbrella sampling can substantially improve sampling efficiency.3,49,50 In practice, however, the difficulty of free-energy reconstruction generally increases with the number of CVs. This dimensionality dependence places strong pressure on identifying a compact representation, often consisting of only one or two coordinates, that retains the essential slow modes of the molecular process of interest.

Identifying such coordinates is itself a challenging problem. From a theoretical perspective, reaction coordinates related to the committor function provide an ideal description of transition progress,51 while recent machine-learning approaches offer increasingly sophisticated strategies for discovering low-dimensional representations from simulation data.32,39,52–57 In both cases, however, determining useful coordinates often requires representative sampling of the relevant metastable states and transition regions. This creates an inherent interdependence between CV discovery and enhanced sampling: effective CVs enable efficient exploration, yet substantial exploration may be required before those CVs can be reliably identified. As a result, the requirement to determine an optimal low-dimensional representation can itself become a bottleneck in free-energy calculations.

An alternative strategy is to reduce the dependence on aggressive dimensionality reduction by developing sampling methods that remain effective in richer CV spaces. Rather than requiring the essential molecular coordinates to be compressed into one or two variables before enhanced sampling begins, such approaches could permit sampling in a larger set of physically meaningful coordinates, with dimensionality reduction or mechanistic interpretation performed subsequently. Existing methods, including metadynamics,20 adaptive biasing approaches,28 and temperature-accelerated molecular dynamics (TAMD),24 can be extended to multidimensional CV spaces, but their efficiency can deteriorate as dimensionality increases. In particular, these trajectory-based enhanced-sampling methods still rely on the progressive exploration of CV space through biased or accelerated dynamics and may remain hindered by slow degrees of freedom or barriers not adequately represented by the chosen coordinates.

Here, we develop an enhanced-sampling framework based on Stein variational gradient descent (SVGD),58 a particle-based variational inference method that transports an ensemble of particles toward a target probability distribution. We adapt SVGD to molecular simulations by defining particles in CV space and coupling each particle to a restrained MD simulation. The restrained simulations provide local thermodynamic information, while SVGD adaptively redistributes the particle locations so that the ensemble approaches the target equilibrium distribution.

This formulation provides a different mechanism for allocating sampling across CV space. Rather than relying on predefined sampling windows or on the progressive exploration of biased or accelerated trajectories, the locations of the restrained simulations are adaptively redistributed through a collective variational transport step. Consequently, exploration does not rely solely on spontaneous barrier-crossing events occurring on their natural timescales. The use of multiple particles also makes the framework compatible with parallel computation, because the restrained simulations associated with different particles can be propagated concurrently between SVGD updates.

We evaluate the method using a series of systems of increasing complexity and CV dimensionality to demonstrate the accuracy and efficiency of CV-SVGD in recovering equilibrium distributions and associated free-energy landscapes. To facilitate adoption of the approach, we additionally provide a Python implementation that couples SVGD updates with restrained molecular simulations and can be adapted to different molecular systems and CV definitions.

Theory

We explore an enhanced sampling strategy based on SVGD, in which an ensemble of particles is deterministically transported toward a target Boltzmann distribution. In the molecular simulation context, the particles can be interpreted as molecular configurations or, as described below, points in a lower-dimensional CV space. The SVGD update combines two effects: a force-like term that drives particles toward regions of high target probability, and a kernel-mediated repulsive term that maintains diversity among particles and prevents collapse into a single mode.

Stein variational gradient descent

SVGD formulates distributional approximation as a deterministic optimization problem. Given a target distribution $p(\mathbf{x})$, the goal is to evolve an initial trial distribution $q(\mathbf{x})$ so that it more closely matches $p(\mathbf{x})$ by minimizing the Kullback–Leibler (KL) divergence.

SVGD considers a smooth perturbation of the current particle positions,

where $\varepsilon$ is a small step size and $\boldsymbol{\phi}$ is a vector-valued velocity field. If $\mathbf{x} \sim q$, then the transformed variable $\mathbf{z} = T(\mathbf{x})$ follows a transformed distribution denoted by $q_T$. The agreement between $q_T$ and the target distribution $p$ is measured by

Taking the derivative of this KL divergence with respect to $\varepsilon$ at $\varepsilon=0$ gives

where $\mathcal{A}_p$ is the Stein operator,

Thus, decreasing the KL divergence is equivalent to choosing a velocity field $\phi$ that maximizes the expected Stein operator trace.

To obtain a tractable closed-form update, SVGD restricts $\phi$ to a norm-bounded vector-valued reproducing kernel Hilbert space (RKHS), $\mathcal{H}^d$. Let $k(x,x')$ be a positive-definite kernel with scalar RKHS $\mathcal{H}$. Each component $\phi_m$ of the vector field $\phi=(\phi_1,\ldots,\phi_d)$ lies in $\mathcal{H}$ and satisfies the reproducing property,

Under this RKHS constraint, the steepest descent direction of the KL divergence is

In practice, the expectation over $q$ is approximated using an ensemble of $n$ particles, $\{\mathbf{x}_i^\ell\}_{i=1}^n$, at iteration $\ell$. The empirical SVGD update is

with

For a Boltzmann target distribution,

Therefore, the first term in Eq. (10) moves particles along a kernel-smoothed thermodynamic force toward regions of lower potential energy and higher Boltzmann probability. The second term, $\nabla_{\mathbf{x}_j} k(\mathbf{x}_j,\mathbf{x}_i)$, acts as a deterministic repulsive interaction among nearby particles, encouraging the ensemble to represent the full target distribution rather than collapsing into a single minimum.

SVGD in collective-variable space

Direct application of SVGD in the full molecular configuration space is challenging for high-dimensional molecular systems. First, the dimensionality of $\mathbf{x}$ can be very large, requiring many particles to represent the target distribution accurately. Second, kernel methods become increasingly difficult in high dimensions because pairwise distances become less informative and kernel bandwidth selection becomes less robust.59

To address these challenges, we apply SVGD in a lower-dimensional CV space. Let

denote a set of $N$ CVs, and let $\boldsymbol{\zeta}$ denote a point in this CV space. The equilibrium probability density in CV space is

where $F(\zeta)$ is the free energy surface and $Z_{\zeta}$ is the corresponding normalization constant. Formally,

where $Z=\int d\mathbf{x}\,\exp[-\beta U(\mathbf{x})]$.

In CV space, the SVGD update has the same form as in configuration space (Eq. (9)), but with $\zeta$ replacing $\mathbf{x}$ and the CV Boltzmann distribution replacing the configurational Boltzmann distribution. The score function becomes

Thus, SVGD in CV space requires estimates of the mean force $\nabla_{\zeta}F(\zeta)$.

We estimate the mean force at a CV-space particle from restrained molecular simulations. For a harmonic restraint centered at $\zeta$,

the corresponding restrained free energy is

Differentiating Eq. (17) with respect to the restraint center gives

where $\langle\cdot\rangle_{\kappa,\zeta}$ denotes an equilibrium average in the restrained ensemble centered at $\zeta$. In the stiff-spring limit, $F_{\kappa}(\zeta)$ approximates the underlying free energy surface up to a smoothing induced by the restraint width. Therefore, the restrained mean force in Eq. (18) provides a practical estimate of $\nabla_{\zeta}F(\zeta)$ for use in the CV-space SVGD update. This estimate was shown by Maragliano and Vanden-Eijnden 24 to have an error of order $O\left(\frac{1}{\beta\kappa}\right)$.

This formulation converts free-energy exploration into an interacting-particle optimization problem in CV space. Each SVGD particle represents a CV location at which local molecular simulations estimate the mean force, while the kernel interaction couples the particles so that they collectively explore and approximate the equilibrium distribution $p(\zeta) \propto\exp[-\beta F(\zeta)]$.

Results

SVGD as an Enhanced Sampling Method

The theory above provides a direct route for applying SVGD to molecular free-energy calculations. In CV-SVGD, each particle represents a location in CV space and is associated with a restrained MD simulation centered at that location. At each iteration, these simulations provide local estimates of the mean force, which are used in the SVGD update to collectively redistribute the particles in CV space. The restraints then guide the underlying molecular configurations toward the updated particle positions, and the procedure is repeated iteratively (Figure 1).

Because practical SVGD calculations can use a relatively small number of particles, the particle positions themselves may provide only a sparse representation of the equilibrium distribution and may therefore be insufficient for directly estimating a detailed free-energy surface. Instead, we reconstruct the free-energy surface from the mean forces evaluated at the particle locations. In this way, SVGD determines where molecular sampling is performed, while the resulting mean-force information is used for free-energy reconstruction.

We implemented CV-SVGD in the Python package CVSVGD. The implementation automates the coupling between restrained molecular simulations, mean-force estimation, and SVGD particle updates, and supports parallel execution of the restrained simulations associated with different particles. The workflow can be applied to arbitrary choices of CVs and iterated for a prescribed number of steps or until the particle distribution has converged.

Validation of CV-SVGD on an Analytically Solvable Model

We first validate CV-SVGD using an analytically solvable three-dimensional double-well model. The potential is defined as

where $\mathbf{x}=(x,y,z)$, $k_1$ and $k_2$ determine the widths of the two wells, and $\mathbf{x}_1$ and $\mathbf{x}_2$ denote their centers. We consider two parameterizations of this model, with different values of $(k_1,k_2)$, that produce relatively low and high free-energy barriers between the two basins. (see Methods).

To mimic the application of CV-SVGD to a molecular system, we treat $(x,y)$ as the CVs and integrate over $z$, thereby defining a two-dimensional CV space. We use 1000 SVGD particles and evaluate the mean forces analytically using Eq. (23).

For the low-barrier system, the SVGD particles are evolved at the physical inverse temperature, $\beta= 1$, and therefore directly target the equilibrium probability distribution. For the high-barrier system, we introduce a simple temperature-scaling strategy to improve coverage of CV space. The mean force $\nabla F(x,y)$ continues to correspond to the free-energy surface of the original system at $\beta= 1$, but its contribution to the SVGD update is scaled by a smaller effective inverse temperature, $\beta_H = 0.024$. The particles therefore evolve toward the flattened distribution

which reduces the contrast between low- and high-free-energy regions and prevents the particle ensemble from becoming overly concentrated within the free-energy minima. This strategy is conceptually analogous to the temperature acceleration used in AFED and TAMD.19,24,27 Importantly, changing $\beta_H$ alters only the target distribution used to guide placement of the SVGD particles; it does not modify the underlying free-energy surface or the mean forces used for its reconstruction.

After convergence, the final SVGD particle distributions closely follow the corresponding analytical target distributions (Figure 2C,F). For the low-barrier system, the particles reproduce the equilibrium distribution at $\beta= 1$. For the high-barrier system, they instead sample the deliberately flattened distribution defined by $\beta_H = 0.024$, again showing close agreement with the analytical result.

We next use the particle locations as centers at which mean-force information is available and reconstruct the free-energy surface by fitting the mean forces with Gaussian radial-basis functions.60 For both the low- and high-barrier systems, the reconstructed surfaces closely reproduce the corresponding analytical free-energy surfaces (Figure 2A,B,D,E).

Together, these results demonstrate that CV-SVGD can use local mean-force information both to adaptively distribute sampling across CV space and to recover the underlying free-energy surface, even when the particle distribution is generated using an accelerated effective temperature.

Aqueous-Phase Alanine Dipeptide

We next apply CV-SVGD to aqueous-phase alanine dipeptide as a proof of concept, using the backbone Ramachandran dihedral angles $\phi$ and $\psi$ to define the two-dimensional CV space (Figure 3A). Alanine dipeptide is a prototypical benchmark for enhanced-sampling methods,60–63 because its conformational free-energy landscape is well characterized and can be obtained accurately using established approaches. We therefore constructed a reference free-energy surface in $(\phi,\psi)$ space using umbrella-sampling simulations.3 Unbiased statistical weights were obtained by solving the multistate Bennett acceptance ratio (MBAR) equations64 using FastMBAR,65 and the resulting reference PMF is shown in Figure 3B.

As in the high-barrier toy model, we used temperature scaling to improve the spatial coverage of the SVGD particles. The restrained MD simulations used to estimate mean forces were performed at the physical temperature of 300 K, while the SVGD particle dynamics used an effective temperature of 5000 K. We propagated 200 particles for 40 SVGD iterations, using 1 ps of restrained MD per particle at each iteration to estimate the local mean force. After the final iteration, the particle locations were fixed and additional 50 ps restrained simulations were performed at each location to obtain more accurate mean-force estimates for free-energy reconstruction.

The final SVGD particle distribution is consistent with the temperature-scaled target distribution derived from the umbrella-sampling reference (Figure 3C). The one-dimensional marginal distributions reproduce the major features of the reference distribution, with residual deviations that likely arise from finite particle number and histogram discretization. More importantly, the free-energy surface reconstructed from the mean-force data closely reproduces the umbrella-sampling result (Figure 3D). CV-SVGD recovers the locations and relative organization of the major metastable basins and the intervening barriers, with an overall root-mean-square error of 1.05 $k_\mathrm{B}T$ relative to the umbrella-sampling PMF. Although the production calculation was propagated for 40 SVGD iterations, convergence analysis indicates that comparable accuracy may already be achieved after approximately 20 iterations (Figure S1).

The complete CV-SVGD calculation required approximately 18 ns of aggregate restrained-MD sampling. By comparison, the reference umbrella-sampling calculation consisted of 625 windows, each simulated for 500 ps, corresponding to 312.5 ns of aggregate MD sampling. Thus, for this benchmark, CV-SVGD reproduced the reference free-energy landscape using substantially less total MD sampling.

Aqueous-Phase Alanine Tripeptide

To further test CV-SVGD in a higher-dimensional setting, we applied the method to aqueous-phase alanine tripeptide using the four backbone Ramachandran dihedral angles ($\phi_1$, $\psi_1$, $\phi_2$, $\psi_2$) as CVs (Figure 4A). Direct construction of a high-resolution four-dimensional reference free-energy surface using conventional umbrella sampling would be computationally demanding. For example, a grid containing 25 windows along each coordinate would require $25^4 = 390{,}625 umbrella windows. We therefore constructed reference data using six independent two-dimensional umbrella-sampling calculations, corresponding to all unique pairwise projections of the four-dimensional CV space.

As for alanine dipeptide, we used temperature scaling to improve coverage of CV space.

The restrained MD simulations used to estimate mean forces were performed at 300 K, whereas the SVGD particle dynamics used an effective temperature of 10,000 K. An ensemble of 200 particles was propagated for 200 SVGD iterations, with 1 ps of restrained MD per particle and per iteration used to estimate the local mean forces. After the final iteration, longer restrained simulations were performed at the converged particle locations to obtain more accurate mean-force estimates for reconstruction of the four-dimensional free-energy surface.

Because a direct four-dimensional umbrella-sampling reference was not available, we assessed the reconstructed CV-SVGD surface through its six unique two-dimensional projections. The projected CV-SVGD free-energy surfaces reproduce the major metastable basins and intervening barriers observed in the corresponding umbrella-sampling references (Figure 4C,D). Some projections exhibit smoother features and modest underestimation of barrier heights. These differences likely reflect a combination of finite particle density, statistical uncertainty in the estimated mean forces, and smoothing introduced during the radial-basis-function reconstruction and subsequent marginalization.

The particle distributions provide an independent assessment of convergence that does not depend on reconstruction of the free-energy surface. As shown in Figure 5A, the one-dimensional marginal distributions of the converged SVGD particles closely reproduce the corresponding temperature-adjusted distributions obtained from the MBAR-reweighted umbrella-sampling data. The Kullback–Leibler divergences between kernel-density estimates of the SVGD and reference distributions decrease rapidly and largely plateau within approximately 100 iterations for all four CVs (Figure 5B). This agreement suggests that part of the residual discrepancy in the reconstructed free-energy surfaces arises from finite mean-force sampling and fitting rather than from failure of the SVGD particle ensemble to approximate the target distribution.

We next examined whether the adaptive placement of particles by SVGD improves free-energy reconstruction relative to a non-adaptive distribution of sampling locations. Figure 4B compares free-energy profiles reconstructed from mean forces evaluated at the converged SVGD particle positions with those obtained using particle positions distributed uniformly throughout the four-dimensional dihedral-angle space. The SVGD-distributed force centers consistently produce better overall agreement with the umbrella-sampling references.

This comparison highlights an important role of the SVGD dynamics beyond distributional sampling: by preferentially allocating particles to thermodynamically relevant regions while maintaining coverage through kernel-mediated repulsion, CV-SVGD provides a more informative set of locations for subsequent mean-force reconstruction.

The alanine-tripeptide calculation also illustrates a distinction between convergence of the SVGD particle distribution and convergence of the reconstructed free-energy surface. Short, 1-ps restrained simulations were sufficient for the SVGD ensemble to approach the reference marginal distributions, whereas substantially longer simulations at the final particle locations were required to obtain mean forces accurate enough for quantitative free-energy reconstruction. Thus, the particle transport itself appears relatively tolerant of noisy local force estimates, whereas accurate reconstruction of a multidimensional free-energy surface places more stringent demands on the precision of the mean forces. Increasing the number of particles provides an additional route to improving the reconstruction: simulations with 500 particles reduce the RMSE of the projected free-energy surfaces to approximately 1.0 $k_{\mathrm{B}}T$ (Figure S2).

In terms of aggregate simulation cost, the 200-particle SVGD propagation requires 40 ns of restrained-MD sampling. Including 200 ns of additional sampling at each final particle location gives a total cost of 0.24 $\mu$s. For comparison, the six two-dimensional umbrella-sampling reference calculations require approximately 0.375 $\mu$s of aggregate simulation time to reach errors below 2.0 $k_{\mathrm{B}}T$ (Figure S3). Importantly, however, these umbrella calculations provide only six two-dimensional projections. A direct four-dimensional umbrella-sampling calculation with 25 windows per dimension would require 390,625 windows; even at 100 ps per window, this corresponds to approximately 39 $\mu$s of aggregate simulation time. CV-SVGD therefore becomes increasingly attractive as the dimensionality of the target CV space grows, because the number of restrained simulations is determined by the number of particles rather than by a tensor-product grid over all CVs.

Conclusions and Discussion

In this work, we introduced CV-SVGD, an enhanced-sampling framework that combines SVGD with restrained MD to explore free-energy landscapes in CV space. By using local mean-force estimates to collectively transport an ensemble of CV-space particles, the method provides an alternative to conventional trajectory-based exploration while retaining an explicit connection to molecular configurations. Across a series of model and molecular systems, CV-SVGD accurately reproduced equilibrium distributions and free-energy surfaces, including multidimensional landscapes.

CV-SVGD shares several conceptual features with TAMD and related approaches.27,33 In both cases, the evolution of CVs is driven by mean forces estimated from restrained molecular simulations, and an elevated effective temperature can be used to facilitate exploration of CV space. The principal distinction is that TAMD propagates auxiliary CVs through a continuous dynamical trajectory, whereas CV-SVGD evolves an interacting ensemble of particles through a collective variational transport step. As a result, exploration in CV-SVGD is not governed solely by sequential barrier-crossing events in CV space. Information from the entire particle ensemble contributes to each update, while the underlying restrained molecular simulations can be carried out concurrently. This combination of collective transport and parallel sampling may become increasingly advantageous as the dimensionality and complexity of the CV space increase.

We anticipate that this capability may be particularly useful for systems in which the relevant slow degrees of freedom are not known a priori. Examples include biomolecular condensates,66–69 chromatin organization,70–74 and other heterogeneous molecular assemblies whose behavior may depend on many coupled structural and compositional coordinates.75 In such systems, identifying a minimal set of optimal CVs can itself be difficult. A sampling framework that remains effective in a richer CV representation could reduce the need to determine this low-dimensional description before enhanced sampling begins. In the present study, we explicitly reconstructed free-energy surfaces containing up to four CVs. This limit was dictated primarily by the difficulty of constructing and validating reliable high-dimensional reference free-energy surfaces rather than by a demonstrated breakdown of CV-SVGD. We plan to determine the practical dimensional limits of the method in future studies.

Several aspects of the present implementation could be improved further. First, we reconstruct the free-energy surface exclusively from mean-force estimates. Accurate mean forces can require substantial sampling at each particle location and may therefore become a significant component of the computational cost. The restrained trajectories, however, contain additional statistical information beyond their mean forces. Combining force-based reconstruction with configuration-based estimators, for example through MBAR-like reweighting or related multistate approaches, could make more complete use of the generated simulation data and potentially reduce the amount of sampling required to achieve a given free-energy accuracy.

Second, the efficiency of CV-SVGD depends on how the particles are distributed across CV space. In the present implementation, we used an elevated effective temperature to flatten the target distribution and prevent excessive concentration of particles within low-free-energy basins. More adaptive strategies may provide better control over particle placement. For example, the free-energy surface accumulated during successive SVGD iterations could be used to construct an evolving bias potential, analogous in spirit to adaptive-biasing approaches.28,40,76 The particle dynamics could then be updated using this progressively refined estimate, allowing sampling effort to be redistributed toward poorly explored regions while retaining the information accumulated in previous iterations. Such an approach could further improve both coverage of high-dimensional CV spaces and the statistical efficiency of free-energy reconstruction.

Methods

General SVGD Protocol

All CV-SVGD calculations were performed using our Python implementation, CVSVGD. For the molecular systems, restrained MD simulations were performed using OpenMM77 to estimate the local mean forces required for the SVGD updates. Because the restrained simulations associated with different particles are independent between SVGD updates, they can be executed concurrently. We used the Parsl Python library78 to manage parallel execution across distributed computing resources.

The SVGD particle positions were updated using the Adam optimizer79 implemented in PyTorch.80 The interactions between particles were evaluated using the inverse multiquadric (IMQ) kernel, which has favorable convergence properties for certain classes of high-dimensional target distributions.81 The kernel is defined as

where the bandwidth parameter $h$ was determined from the median pairwise particle distance at each iteration. The Adam learning rate was set to 0.6 for all calculations, with additional optimizer parameters specified below.

Following convergence of the SVGD particle ensemble, additional restrained simulations were performed at the final particle locations to obtain more accurate mean-force estimates for free-energy reconstruction. The resulting mean forces were fitted using radial-basis functions (RBFs), following the approach of Maragliano et al. 60``. Details of the RBF fitting procedure are provided in the Supporting Information.

Two-dimensional Toy System

We validated CV-SVGD using an analytically solvable three-dimensional double-well potential, defined in Eq. (19). Taking $(x,y)$ as the CVs and integrating out the $z$ coordinate, the corresponding free-energy gradient is

At $\beta=1$, this expression can be evaluated analytically:

where $\mathbf{s}=(x,y)$ and $\mathbf{s}_1=(x_1,y_1)$ and $\mathbf{s}_2=(x_2,y_2)$ denote the projected well centers.

We considered two parameterizations with different barrier heights. The low-barrier system used $k_1=k_2=2$, whereas the high-barrier system used $k_1=50$ and $k_2=25$. For both systems, the well centers were fixed at $\mathbf{x}_1=(1,1,1)$ and $\mathbf{x}_2=(-1,-1,-1)$.

At $\beta= 1$, the analytical equilibrium distribution in the two-dimensional CV space is

For the low-barrier system, SVGD was performed at $\beta= 1$, targeting the equilibrium distribution $P(x,y)$. For the high-barrier system, we employed temperature scaling with an effective inverse temperature $\beta_H = 0.024$, while retaining the free-energy gradients evaluated at $\beta= 1$. The corresponding SVGD target distribution is

where $\beta_H$ controls the degree of flattening of the target distribution without modifying the underlying free-energy surface.

For both systems, 1000 SVGD particles were initialized uniformly over the region $x,y \in[-5,5]$ and propagated for 150 iterations using the analytical free-energy gradients. The particle positions obtained after the final iteration were used as centers for RBF fitting to reconstruct the free-energy surfaces.

Alanine Dipeptide

CV-SVGD Simulations

For aqueous-phase alanine dipeptide, we initialized 200 SVGD particles uniformly over the two-dimensional Ramachandran dihedral-angle space $(\phi,\psi)$. Initial molecular configurations were generated using tleap in AmberTools.82 Before initiating the SVGD updates, each configuration was subjected to a 10 ps restrained MD simulation to guide the dihedral angles toward their assigned particle positions.

All MD simulations were performed using OpenMM with the ff14SB protein force field83 and the TIP3P-FB explicit water model.84,85 The system was solvated with 423 water molecules in a cubic periodic box of side length 2.4 nm. Long-range electrostatic interactions were treated using the particle mesh Ewald method,86 with a nonbonded cutoff of 1.0 nm. Simulations were propagated using the Langevin middle integrator at 300 K, with a friction coefficient of 5 ps$^{-1}$ and a time step of 1 fs. Molecular configurations were recorded every 0.01 ps for mean-force estimation.

To constrain molecular configurations near the corresponding SVGD particle positions in CV space, we applied a harmonic restraining potential of the form

where $\mathbf{x}$ denotes the molecular coordinates, $\phi(\mathbf{x})$ and $\psi(\mathbf{x})$ are the instantaneous backbone dihedral angles, and $(\phi_i,\psi_i)$ specifies the position of the $i$th SVGD particle in CV space. The restraint force constant was set to $\kappa= 2500$ kJ mol$^{-1}$ rad$^{-2}$.

The SVGD particle ensemble was optimized at an effective temperature of 5000 K, while the restrained MD simulations used for mean-force estimation remained at 300 K. Particle positions were updated using the Adam optimizer with $\beta_1 = 0.6$, $\beta_2 = 0.999$, $\epsilon= 0.1$, and a learning rate of 0.6.

Following each SVGD update, an additional 0.1 ps restrained MD simulation was performed to allow the molecular configurations to relax toward the updated particle positions. This intermediate relaxation step reduced the occurrence of excessively large forces arising from abrupt changes in restraint centers, allowing the use of larger SVGD updates without compromising numerical stability.

After the final SVGD iteration, the particle positions were fixed and an additional 50 ps restrained MD simulation was performed at each location to obtain more accurate mean-force estimates for free-energy reconstruction. The free-energy surface was reconstructed by RBF fitting using von Mises functions, which account for the periodicity of the dihedral angles. The one-dimensional von Mises basis function takes the form

where $\lambda$ controls the width of the basis function. These periodic basis functions were used in place of Gaussian functions for fitting the free-energy surface in Ramachandran space.

Umbrella Sampling Simulations

Reference free-energy surfaces were obtained using two-dimensional umbrella sampling in the $(\phi,\psi)$ dihedral-angle space. The molecular model and simulation parameters were identical to those used for CV-SVGD, except that the restraint force constant was reduced to $\kappa= 500\ \mathrm{kJ\,mol^{-1}\,rad^{-2}}$ and the integration time step was increased to 2 fs.

Each dihedral angle was discretized into 25 equally spaced umbrella windows spanning $[-\pi,\pi)$, yielding a total of $25 \times25 = 625$ windows. Each window was simulated for 0.5 ns, with configurations recorded every 1 ps.

Alanine Tripeptide

CV-SVGD Simulations

For aqueous-phase alanine tripeptide, we initialized 200 SVGD particles uniformly over the four-dimensional backbone dihedral-angle space $(\phi_1,\psi_1,\phi_2,\psi_2)$. The SVGD particle ensemble was propagated for 200 iterations using an effective temperature of 10,000 K, while the restrained MD simulations used for mean-force estimation were performed at 300 K. Each restrained simulation was 1 ps in duration.

The molecular system was solvated with 444 water molecules in a cubic periodic box of side length 2.5 nm. All other molecular simulation parameters, including the force field, water model, integrator, and restraint protocol, were identical to those used for alanine dipeptide.

Following the final SVGD iteration, additional 2 ns restrained MD simulations were performed at each particle location to reduce statistical noise in the mean-force estimates. These refined mean forces were subsequently used to reconstruct the four-dimensional free-energy surface through radial-basis-function fitting with von Mises kernels.

Umbrella Sampling Simulations

To validate the four-dimensional free-energy surface obtained from CV-SVGD, we performed six independent sets of two-dimensional umbrella-sampling simulations, corresponding to all unique pairs of the four backbone dihedral angles.

The molecular model and simulation parameters were identical to those used for CV-SVGD, except that the restraint force constant was reduced to $\kappa= 250\ \mathrm{kJ}\ \mathrm{mol}^{-1}\ \mathrm{rad}^{-2}$ and the integration time step was increased to 2 fs. Each dihedral angle was discretized into 25 equally spaced umbrella windows spanning $[-\pi, \pi)$, yielding $25 \times25 = 625$ windows per CV pair and 3,750 windows across all six umbrella-sampling calculations. Each window was simulated for 0.5 ns, with configurations recorded every 1 ps.

Acknowledgement

This work was supported by the National Institutes of Health (Grant R35GM1335w80). I.R. was the recipient of the National Science Foundation Graduate Research Fellowship.

Competing interests: The authors declare that they have no competing interests.

Data and materials availability: CVSVGD Python API is available at: https://github.com/ZhangGroup-MITChemistry/CVSVGD

Supporting Information: Technical Details of our RBF Fitting process, convergence plots for the alanine dipeptide system, repeat of figure 4 with 500 CV-space particles, convergence plots for umbrella sampling simulations, as well as videos of CV-space trajectories for the two toy model systems, the alanine dipeptide system and alanine tripeptide system.

Supporting Information

Fitting Radial Basis Functions to Mean-Forces

We utilized the radial basis function (RBF) fitting method of Maragliano et al.S1 to recover the free energy surfaces for the systems investigated. The fitting method is described in the following. An RBF is of the following form:

where $k(\cdot)$ is any suitable basis function, $\zeta= (\xi_1, \ldots, \xi_d)$ is a vector of $d$ CV-space coordinates, $\alpha_i$ is a constant which assigns a weight to the $i$-th RBF component, $\mu_i$ is the center of the $i$-th RBF component, and $\sigma$ is a bandwidth parameter controlling the basis function’s width. In our investigations we use the typical gaussian kernel for non-periodic CVs:

and the von Mises kernel for periodic ones:

Given a set of mean forces $\{\nabla_{\zeta}F_i\}_{i=1}^{N}$, and corresponding force centers $\{\zeta_i\}_{i=1}^{N}$, we wish to parameterize our RBF $\phi(\zeta)$ such that the following loss is minimized:

The RBF’s gradient $\nabla_{\zeta}\phi(x)$ is a sum of the gradient of basis functions expanding as follows:

The loss function can be then be equivalently rewritten as a matrix equation:

Where $K$ is a $(dN) \times N$ matrix, where elements of matrix $K$ are given by

$\alpha= (\alpha_1, \ldots, \alpha_N)$ is a vector of the weight parameters, and

is a vector containing all mean-forces component-wise, in other words a flattened $d \times N$ matrix of the combined mean-forces. To prevent over-fitting, we also employ regularization on $\alpha$ with a regularization parameter $\lambda> 0$, giving the final loss of:

A $\lambda= 10^{-9}$ was used for all RBF fittings. The optimal solution to eq. (33) can be found using reduced QR-factorization, where we factorize the augmented matrix $A$:

Where $I_N$ is an $N \times N$ identity matrix. Factoring $A = QR$, gives:

Where:

and $\vec{0}$ is a N-component $0$-vector. Using eq. (34), the optimal $\alpha$ can be found via back-substitution. To ensure the optimal fit of all parameters $\alpha$ and $\sigma$, we use a PytorchS2 implementation of the described fitting procedure, allowing the $\sigma$ to be additionally optimized via gradient descent. This is done by initializing the fitting procedure with a guess $\sigma$ value, finding the optimal $\alpha$ given that $\sigma$ value, and then computing the gradient of eq. (33) with respect to $\sigma$, and iterating until convergence. Notably, reduced QR-factorization is fully differentiable in the Pytorch implementation, enabling quick convergence to the optimal parameters. Additionally, it allows for anisotropic kernels via a unique $\sigma$ per CV dimension to better fit the free energy surface.

Supplemental Figures