Figure compares the accuracy of the semi-analytic and numerical methods with a fixed simulation run length. The mean and standard error for each gradient is computed from five independent runs using the TIP3P parameters with a simulation length of 20ns. Finite difference step sizes of δk λ = 0.001, 0.01 and 0.1 in the mathematical parameters were tested for both methods. When numeric gradients were used, the statistical errors were largest for step sizes of 0.001 and smallest for 0.01. For σ LJ , increasing the step size to 0.1 resulted in a different mean and larger standard error, indicating this step size was outside the linear regime; we did not observe this for q H and LJ . The semi-analytic gradients are computationally less costly but also have higher uncertainty than the numerical gradients, thus we concluded that numerical gradients with a step size of 0.01 provide the most statistically precise surface tension gradients for a fixed simulation length. These conclusions are based on our choice of the prior widths for the parameters; for a different choice of prior width, the recommended step size may be obtained by ensuring the step size in physical parameters, i.e. δK λ = δk λ p λ matches the current results. Taking the total computational cost into account, we compared the numeric gradients from 5 ns simulation with the semi-analytic gradients from 20 ns simulation, both with the optimal step size 0.01, as shown in Figure . Two sets of parameters were used, namely the TIP3P parameters and the TIP3P-ST parameters. The results show that for the TIP3P parameters, the numeric and semi-analytic gradients agree relatively well with comparable standard errors. However, when evaluated at the final TIP3P-ST parameters, the semi-analytic gradients have larger errors than the numeric gradients for the q H and σ LJ parameters, while LJ exhibits the opposite behavior. The small errors for the semi-analytic gradients of LJ may be due to the intrinsically small value of the gradient (i.e. in the limit of infinite simulation time). An intrinsically small gradient would reduce the error bars of the analytic gradient but not the numerical gradient, as the latter contains statistical noise from independent estimations of the surface tension and contributes a constant term to the error. The scale-independent behavior of the numerical gradient error is confirmed by comparing the standard error across parameters; for TIP3P these errors are (15.2, 18.0, 13.3) for q H , σ LJ , LJ respectively, and for TIP3P-ST the errors are (59.0, 32.9, 41.4). The standard error for surface tension gradients are larger overall for TIP3P-ST compared to TIP3P, which may be due to the slower dynamics of the model causing slower convergence of the property. Based on our observation that the statistical errors were mostly smaller using numeric gradients, we decided to use numeric gradients for optimizing the water models in this study.