Regularization was used in the optimization to achieve smooth convergence and to prevent the parameters from moving too far from the starting point to potentially unphysical local minima. For this optimization, the starting point was the 1.3.0 force field. Regularization with a harmonic term can be seen in a Bayesian sense as imposing a Gaussian prior on each parameter with mean of the starting point, and standard deviation equal to the regularization scale 𝜎 𝑝 . We used a data-driven approach to determine the regularization scales used in the fitting procedure. The distribution of parameters in SMIRNOFF99Frosst for each parameter type was first plotted. Since the distributions were not bell-shaped (as seen in SI Figure . The overall optimization of the force field through this process was considered to have converged when it satisfieed at least two out of the following three criteria: total objective function value, including the regularization penalty, to reach a value ≤ 0.1; the norm of gradient on parameters to reach a value ≤ 0.1; and the optimization step size to be ≤ 0.01.