The procedure for evaluating surface tension was implemented into the ForceBalance automated parameter optimization software, which uses the OpenMM library to carry out the NVT and NPT MD simulations, thus allowing the entire optimization procedure to be carried out in a single reproducible calculation. Although thermodynamic fluctuation formulas were used to estimate the parametric derivatives of thermodynamic properties simulated in the NPT ensemble in previous parameterization of TIP3P-FB and TIP4P-FB models, we found that in the case of surface tension, the parametric derivatives estimated in this way contained such high levels of statistical noise that it was more efficient to calculate parametric derivatives numerically via a 3-point finite difference formula, which involved running two separate simulations for each parameter being optimized. Details of the error analysis are described in Section 3.3.