As discussed before, this large-scale parameterization of water models is made possible by the previous two approaches which dramatically increase the efficiency of free energy calculations for a very large number of thermodynamic states with little or no configurational overlap by using explicit maps. Like all other parameterizations, our scheme also involves iterations to convergence. instead of gradient methods but such more robust direct methods converge slowly. We cannot use simplex or gradient-based method because the simplex algorithm does not support multi-objective optimization and in gradient-based methods the use of the gradient with very small step length, when optimal step length is not known, makes the convergence very slow. Calculation of optimal step length or the use of Hessian-based methods is computationally prohibitive, though in theory reweighting could be used to construct the gradient. Instead, we use a mesh adaptive direct search method which is well-suited for our problem type as it avoids the bottlenecks associated with the gradient based methods. We discuss the efficacy of different optimization methods and compare them further in the methods section. As discussed above, we replace the time-limiting step of MD simulations by a faster reweighting plus configuration mapping step to estimate properties at each iteration. The configurational mapping reduces the number of samples required, which must be generated using the initial set of parameters, and the reweighting reduces the need to generate new samples at new iterated sets of parameters.