For the solvated ligand and protein-ligand complex systems, a λ-schedule of 9 equally distributed windows was used (0.1, 0.2, 0.3, …, 0.9). Gradient means and variances were calculated using the bootstrap method. For the addition of the virtual bond restraints, 7 unequally distributed λwindows were used (0.0, 0.05, 0.1, 0.2, 0.3, 0.5, 1.0). Free energies for the ligand, the protein, and the restraint addition were obtained via the trapezoid rule. The free energy of adding virtual bond restraints for the non-interacting ligand was calculated using the Boresch formula. A total of 39 independent simulations were performed for all protocols. All uncertainties were quantified using samples of multiple replicates.