We initialized the NEQ simulations with 500 conformations randomly drawn (with replacement) from 3 independent 5 ns equilibrium simulations (3×4, 000 samples) with the reference and target potential. The switching protocol consists of T perturbation and propagation steps performed alternately. The perturbation step alters the system according to the con-trol parameter λ, and the propagating step alters the coordinates at fixed λ. For a fixed annealing length T , the potential is scaled with λ = t/T at switching time t between reference U MM and target U NNP potential, such that U (t) = (1 -λ)U MM + λU NNP . The propagation step consists of a single integration step at fixed U (λ) using Langevin dynamics with a 1 fs timestep. Dimensionless work values (w) were accumulated throughout the annealing protocol. We tested different annealing lengths T (5, 10, 20, 50 ps) and used the EXP (exponential averaging) estimator to recover the free energy estimates for unidirectional work distributions and the BAR (Bennett's acceptance ratio) estimator to obtain the free energy estimate from bidirectional work distributions as implemented in the pymbar package. Single-step free-energy perturbation FEP was performed using 12,000 samples from the accumulated 3 × 5 ns equilibrium simulations at λ = 0 or λ = 1. The potential energy difference ∆U is calculated using