Free energy is a thermodynamic observable of importance for protein-ligand complexes in the drug discovery context. Therefore, we also look at the extent of variation in free energy estimates obtained using independent "long" replicas of MD simulations. We used ∆G direct as a measure of absolute binding free energy which was originally developed by De Jong et al., and later applied to 10 microsecond trajectories by Pan et al. in 2017 (details in section 3). In the present work, we demonstrate that ∆G direct varies substantially between separate independent long timescale replicas, and hence once again individual "long" simulations do Figure : Running averages of ∆G direct for tofacitinib binding to ADRP at all of the identified ADRP binding sites (top two panels) and for the experimental binding sites of 3CLPro systems (bottom two panels) for all ten or twelve 10µs trajectories. The horizontal black dashed lines in the plots of 3CLPro systems correspond to the respective experimental binding affinities. not provide reliable binding free energy estimates. The salient point here is that, contrary to received wisdom in the literature on molecular dynamics, averaging over an individual long timescale simulation is not equivalent to averaging over an ensemble of simulations. Indeed, thermodynamic quantities arise from ensemble averaging in statistical mechanics and unless one averages over a timescale of the order of a Poincaré recurrence, a one-off MD trajectory will produce the wrong results. Compounding this, a one-off simulation does not provide the means to compute precise results or conduct meaningful uncertainty quantification. Figure shows the running averages of ∆G direct for all four binding sites of the ADRP system (top two rows) as well as for crystallographic sites of all 3CLPro systems (bottom two rows) from all replicas that sample them. The inter-replica variation is clearly visible from these plots for all systems. This variability shows that results obtained from individual "long" trajectories are not reproducible or precise. In fact, it is evident that ∆G direct estimator does not even converge for some replicas, even though it does for others. This itself is a source of uncertainty and a strong motivation to perform ensemble simulations. Nevertheless, even when all or most replicas do converge, their predicitions vary non-negligibly. Figures displaying running averages of ∆G direct for all binding sites of all systems studied have been included in the Supporting Information (Figure ). All of them show behaviour similar to that discussed above in terms of ∆G variability. Table includes the mean ∆G direct values along with error bars for the most frequently visited binding site (which is not the experimental binding site for PLPro-GRL and 3CLPro-RQN (with LZE)) across all replicas for each system. It also provides the spread (that is, the difference between the two extreme values) for each such binding site which is around 2-3 kcal/mol for most cases but can be as high as 7 kcal/mol (for instance 3CLPro-93J). Another point worth noting from Figures and is that the first binding event occurs at varying time durations across ensemble members as captured by the different onset simulation times of the running average plots.