. There is no free energy cost to remove the ligand from the binding pocket in the last step in Fig. (5 → 6). Software. All MD simulations were performed using the GPU-accelerated pmemd.cuda module of AMBER20. All 𝑑𝑉 𝑑𝜆 gradient time series data (Eq. 3) were extracted with the alchemlyb python package. Decorrelation was performed using the pymbar python package time series analysis module, whereby decorrelated samples were obtained by subsampling with the statistical inefficiency rounded up to the nearest integer value. All hydrogen bonding interactions between the protein and ligand were analyzed with CPPTRAJ. All input coordinates, topology and parameters for conventional MD simulations were obtained using Ambertools18. Molecular Systems Simulated RBFE CDK2 system setup and TI simulations. Simulation starting structures for protein-ligand complexes were extracted from the GitHub repository . For both the protein-ligand complex and solvated ligand mutations, two separate λ-schedules were used. The first employed the following twelve λ-windows: steps of minimization with the gradient descent method; 2) 50 ps of heating from 0.1 K to 300 K in the NVT ensemble; 3) 300 ps of density equilibration in NPT; and 4) production simulations in NVT with gradient averages obtained via the bootstrap method. Harmonic RMSD restraints were imposed on heavy atoms of the protein and the ligand during minimization and heating and were gradually removed during density equilibration; no restraints were used during production simulations. Special care was given to the λ=0.98408 window of the nine-point λ-schedule to avoid numerical instability. For this window, the structure obtained after the first 1000 steps of minimization with gradient descent of the λ=0.91802 window was used as the input for the second 1000 steps of minimization, and then the protocol proceeded as normal. Soft-core potentials were utilized for all simulations. For the 12-window λ-schedule production simulations, which was implemented to closely mimic the production protocol of Song et al., a 1 fs timestep was used without SHAKE and a Berendsen thermostat was employed. All other protocols utilized a Langevin thermostat. For the 9-point λ-schedule, which was implemented according to our standard practices, a 2 fs timestep was used with SHAKE. Free energies of both the protein-ligand complex and solvated ligand alchemical steps were obtained using the gaussian quadrature rule.