Classical force field parameterization is well-known to be computationally intensive. Traditionally, bonded parameters such as bond length, bond angle and dihedral angle, are estimated using quantum chemical calculations and/or fit to spectroscopic data. However, the nonbonded parameters, especially Lennard-Jones (LJ) parameters, are usually difficult to determine from first principles calculation or from experiments. The nonbonded parameters, charge q and LJ σ, are instead generally fit to reproduce experimental observables. Fitting to experimental thermodynamic observables involves carrying out sufficiently long molecular simulations to estimate observables at each iteration, which require significant computational time. New sets of parameters are generally either guessed based on a trial and error approach or more occasionally by simplex or gradient-based optimization, though other approaches such as optimization using surrogate models have been tried. The time spent in all of these approaches is dominated by the evaluation of the objective function, which requires a fresh generation of samples for each new guess of the parameter set. We can accelerate the force field parameterization process if we can significantly reduce the time spent in sampling at each iteration.