Estimating Free-Energy Surfaces and Their Convergence from Multiple, Independent Static and History-Dependent Biased Molecular-Dynamics Simulations with Mean Force Integration

Addressing the sampling problem is central to obtaining quantitative insight from molecular dynamics simulations. Adaptive biased sampling methods, such as metadynamics, tackle this issue by perturbing the Hamiltonian of a system with a history-dependent bias potential, enhancing the exploration of the ensemble of configurations and estimating the corresponding free energy surface (FES). Nevertheless, efficiently assessing and systematically improving their convergence remains an open problem. Here, building on mean force integration (MFI), we develop and test a metric for estimating the convergence of FESs obtained by combining asynchronous, independent simulations subject to diverse biasing protocols, including static biases, different variants of metadynamics, and various combinations of static and history-dependent biases. The developed metric and the ability to combine independent simulations granted by MFI enable us to devise strategies to systematically improve the quality of FES estimates. We demonstrate our approach by computing FES and assessing the convergence of a range of systems of increasing complexity, including one- and two-dimensional analytical FESs, alanine dipeptide, a Lennard–Jones supersaturated vapor undergoing liquid droplet nucleation, and the model of a colloidal system crystallizing via a two-step mechanism. The methods presented here can be generally applied to biased simulations and are implemented in pyMFI, a publicly accessible, open-source Python library.

In this section, the pyMFI library [1] is displayed, and interested readers are encouraged to run the code themselves.First, the one-dimensional analytical surface (illustrated in Figure 1a) is simulated with Langevin dynamics, and the results are analysed with MFI.Next, the same one-dimensional surface is simulated and analysed using different simulation parameters and the results from the two simulations are patched.After that, a two-dimensional analytical surface (illustrated in Figure 2) is simulated and analysed with MFI twice, the first time employing only a metadynamics bias and the second time biasing the simulation with metadynamics and harmonic potential.Later we report the simulation setup for the one-and two-dimensional system, followed by the setup for the simulations of alanine dipeptide, the nucleation of the supersaturated Argon vapour, and the nucleation of the colloidal system.The Python code is kept at a level accessible to most readers, but interested readers can have a deeper look at the full library at https://github.com/mme-ucl/MFIand make use of the documentation for a better understanding of the functions.Additional Jupyter notebooks available on GitHub provide additional examples.Langevin simulations can easily be performed on a personal computer with the PLUMED software [2].The simulation of molecular systems, such as alanine dipeptide, requires additional software, such as GROMACS [3], to be installed and patched with plumed.However, existing simulations can also be analysed using pyMFI, requiring only the trajectory data represented as CVs and the biasing history (e.g. the COLVAR and HILLS file).

A. 1D Analytical Surface
The first simulation will be performed on a one-dimensional analytical potential presented in figure 1  .It will be simulated with a Langevin dynamics simulation with a Metadynamics bias.The trajectory data represented as CVs (name of file: "position") and the biasing history (name of file: "HILLS") is analysed with MFI to determine the force-terms (i.e.average mean force, probability density, standard error of the mean force, and other relevant results).Lastly, the average mean force is integrated to obtain the FES, and the results are shown in figure S1. # Compute the time -independent mean force results = MFI1D .MFI_1D ( HILLS = HILLS , position = position , bw =0.03 , kT =1 , log_pace =4000 , error_pace =200 , min_grid = -2 , max_grid =2 , nbins =401 , WellTempered =1 , u se _ w ei g h ted_st_dev = False ) X , Ftot_den , Ftot_den2 , Ftot , ofv_num , FES , ofv , ofe , cutoff , error_evol , f e s _ e r r o r _cu toff _evol = results # integration on a non -periodic domain

B. Combining Two Simulations
Having performed and analysed the first simulation, the force terms are saved so that they can later be combined with the next simulations.Next, an additional simulation is conducted on the same system, this time changing some parameters such that the simulation starts in the basin on the right and the metadynamics of Gaussian hills are higher.After the simulation, the trajectory and biasing history are analysed with MFI, combined with the previous force terms, and the Mean Force is integrated to find the FES.The results are shown in figure S2. # integration on a non -periodic domain The next simulation will be performed on the two-dimensional double well analytical potential presented in figure 3 a), which is defined as F exact (s 1 , s 2 ) = 1.35sHaving analysed the simulation on the two-dimensional double well, one can see in figure S3 b) that the transition region is not sampled very well, and the error of the mean force, depicted in figure S3 c), is relatively high in the transition region.Consequently, one might want to sample the transition region more extensively.In that case, a harmonic potential is applied in the centre (x, y = −0.5, −0.5, with κ x = κ y = 40kJ/mol) while also employing MetaD, to provide a consistent sampling of that region.Similar to the one-dimensional example, the second simulation is analysed, and the force terms of both simulations are combined.It is possible to use the optional argument "base terms" in the MFI 2D function, which requires the force terms of a previous simulation as input.If this option is used, the on-the-fly error calculated is directly patched with the force terms specified in "base terms", yielding the combined error.Also, the ofe history of both simulations can be merged, resulting in a continuous error progression.However, the force terms returned by the MFI 2D function result from the "HILLS" and "position" data and need to be patched with the other force terms using the patch 2D function.

S1. CONVERGENCE ESTIMATOR
The convergence estimator introduced in section III.measures the statistical variance (or standard error) of the mean force.By following its progression, one can better understand which regions of the CV-space have converged, and which regions require further sampling.However, the calculation of the mean force (see Eq. 2 and Eq. 3) depends on the probability density (see Eq. 4), which in turn depends on h, the bandwidth of the Gaussian kernel.Therefore, the choice of the bandwidth affects the standard error of the mean force.Choosing a bandwidth that is too small results in a noisy probability density, leading to a noisy estimation of the mean force and the free energy surface.Consequently, the standard error of the mean force tends to be higher with lower bandwidths.Conversely, a bandwidth that is too large, will result in a probability density that is excessively smooth, impairing its ability to capture local changes in sampling accurately.However, this smoothing leads to a smoother average mean force, and thus, a lower standard error of the mean force.Nonetheless, the convergence estimation does not provide an accurate error, but rather a qualitative assessment of the estimate of the average mean force.A lower bandwidth will result in a higher standard error of the mean force and vice versa, but the qualitative trend in the progression of the standard error remains consistent, as shown in Figure S5.Ultimately, the optimal choice for the bandwidth depends on the roughness of the free energy surface, and the amount of sampling available for the estimation of the probability density.

A. Alanine Dipeptide
The conformational structures of alanine dipeptide in vacuo were obtained from molecular dynamics simulations under the effect of a well-tempered metadynamics bias, with the two backbone torsion angles Ψ and Φ as collective variables.The molecular dynamics simulations were performed with GROMACS [3] 2021 patched with PLUMED[2] 2.7.5, using the AMBER99SB force field [4].The system consists of one alanine dipeptide molecule, simulated in a box of size 2.86 × 2.86 × 2.02 nm with no periodic boundaries.Prior to the production run, energy minimization was performed with a steepest descent algorithm with a residual force tolerance of 10 kJ × mol −1 nm −1 .A cut-off of 1.2 nm for non-bonded interactions was chosen, long-range intermolecular interactions are accounted for using the Particle Mesh Ewald (PME) approach and all bonds are constrained using the LINCS algorithm.The molecular dynamics simulation was executed under NVT conditions (canonical ensemble) at 300K, using a modified Berendsen thermostat (V-rescale thermostat) with a leapfrog integrator and a time step of 2 fs for a total simulation time of 40 ns.The well-tempered metadynamics bias was constructed by adding a Gaussian hill every 2 ps, with a height of 0.5 kJ/mol, a bandwidth of 0.2 for both Ψ and Φ, and a factor of 5.
Another simulation was performed under the same conditions for 500 ns.This simulation time was long enough that the resulting FES was fully converged and used as a reference surface for the calculation of the average absolute deviation.
The input files for these simulations are available on PLUMED-NEST.The forward and backward transitions of a supersaturated Argon vapour to a liquid droplet were simulated multiple times with GROMACS 4.6.3patched with PLUMED[2] 2.0.The system consists of 512 argon atoms, and their interactions were described by a Lennard-Jones potential with ϵ = 0.99797 kJ/mol and σ = 0.3405 nm.The Lennard-Jones potential was truncated at the cutoff length of 5 σ.The molecular dynamics simulation was executed under NVT conditions (canonical ensemble) at a constant temperature of 72K, controlled using the Bussi-Donadio-Parrinello thermostat [5].The time step for the integration of the equations of motion was set to 5 fs, and the simulation was allowed to run until the energy barrier was crossed and the other metastable state was sampled.This setup was conducted 50 times for the forward transition (condensation) and 50 times for the backward transition (evaporation) and repeated for supersaturation levels of 11.95, 14.03, 15.57 and 16.86, which was set by specifying the box dimension l, shown in table S1.Table S1 also includes the initial pressure p, the width of the deposited Gaussians δ, their initial height ω 0 , the deposition stride of the Gaussians ∆t, and the biasfactor γ.For the lowest level of supersaturation (S 1 = 11.95), a large and steep energy barrier was expected so that an additional wall potential was used for the forward transition.

C. Two-step Colloidal Model Nucleation
Molecular dynamics simulations of the colloidal system were performed in the canonical ensemble, tempered to 2T * with 421 particles in a cubic box of length 92.83σ using LAMMPS [6].The colloidal particles are modelled via a Derjaguin-Landau-Verwey-Overbeek (DLVO) potential [7][8][9] with a cutoff of 12.5σ.Details of the potential, its expected thermodynamic behaviour, and additional simulation details are available to reference [10].
The n and n(Q6) CVs are used to describe the nucleation mechanisms of the colloidal system, which were approximated with a graph neural network.Out of the four simulations, three were biased with WTmetaD, each employing a different bias factor (γ 1 =40, γ 1 =50, γ 1 =60), while the fourth simulation was biased with non-tempered metaD.The PLUMED input files are available on PLUMED-NEST (https://www.plumed-nest.org/,plumID:23.026).All machine learning models were trained using the NNucleate package, which is built on top of PyTorch [11].It utilizes functionalities from the MDTraj [12] and MDAnalysis [13,14] packages to train CVs, augment and manage datasets, analyze models, and translate models into PLUMED-readable CV files.These Python scripts are used as CVs through the PLUMED2 fork PyCV [15].Converting the models into a format supported by PyCV involves Alphabet's Jax and Flax packages, and the necessary gradients are obtained using Jax's autodifferentiation implementation [16,17].
a), which is defined as F exact (s) = s 8 − 103e − (s+1.5 FIG. S1. a) MFI estimate of a monodimensional multi-well free energy surface obtained postprocessing a single metadynamics simulation (orange line) with the exact FES (blue line).b) MFI estimate of the probability density of single metadynamics simulation.

#
FIG. S2. a) MFI estimate of a monodimensional multi-well free energy surface obtained post-processing two metadynamics simulations (orange line) with the exact FES (blue line).b) MFI estimate of the probability density of first metadynamics simulation (blue line), second metadynamics simulation (orange line), and combined metadynamics simulations (green line)

#
FIG. S3. Results obtained post-processing a Langevin Dynamic simulation of a two-dimensional double-well system biased with MetaD.a) MFI estimate of the FES.b) Biased probability density.c) Error of the mean force.d) Progression of the error of the mean force (blue line) and progression of the average absolute deviation to the analytical surface (red line)

#
FIG. S4. Results obtained post-processing two Langevin Dynamic simulation of a two-dimensional double-well system.The former simulation was biased with MetaD (results shown in figure S3 and the latter biased with MetaD and US.a) MFI estimate of the FES.b) Biased probability density.c) Error of the mean force.d) Progression of the error of the mean force (blue line) and progression of the average absolute deviation to the analytical surface (red line)

FIG. S5 .
FIG. S5.Comparing the progression of the standard error of the mean force σE divided by the sampled volume v(x-axis) with the average absolute deviation divided by the sampled volume (y-axis).The error progressions in panel a) are calculated with a bandwidth, h, of 0.015.In panel b) a bandwidth of 0.03 was used and in panel c) a bandwidth of 0.06.The blue-red color shift indicates simulation progression.