Computing the free energy of protein-ligand binding by employing molecular dynamics (MD) simulations is becoming a valuable tool in the early stages of drug discovery. However, the cost and complexity of such simulations are often prohibitive for high-throughput studies. We present an automated workflow for the thermodynamic integration scheme with the "on-the-fly" optimization of computational resource allocation for each λ-window of both relative and absolute binding free energy simulations. This iterative workflow utilizes automatic equilibration detection and convergence testing via the Jensen-Shannon distance to determine optimal simulation stopping points in an entirely data-driven manner. It is broadly applicable to multiple free energy calculations, such as ligand binding, amino acid mutations, and others, while utilizing different estimators, e.g., free energy perturbation, BAR, MBAR, etc. We benchmark our workflow on the well-characterized systems, namely, cyclin-dependent kinase 2 and T4 lysozyme L99A/M102Q mutant, and the more flexible SARS-CoV-2 papain-like protease. We demonstrate that this proposed protocol can achieve more than 85% reduction in computational expense while maintaining similar levels of accuracy compared to other benchmarking protocols. We examine the performance of this protocol on both small and large molecular transformations. The cost-accuracy tradeoff of repeated runs is also investigated.