It has been reported in numerous published studies that MD based free energies are sensitive to their initial conditions. Thus, we get a distribution of free energies on performing ensembles of MD simulations with identical input except their starting velocities due to the aleatoric uncertainty in MD. We have shown that aleatoric uncertainty dominates all other forms of uncertainty, and that ensemble simulations are essential. Parametric uncertainty is dampened between input and output, further exacerbating the role of chaos. An additional source of systematic error for chaotic systems originates from the use of floating point numbers. This is true irrespective of the free energy method used including the alchemical ones. Figure displays the distributions of both potential energy derivatives at an arbitrary λ window as well as ∆∆G values for a couple of ligand pairs bound to MCL1 using 40 replicas per λ window. It can be seen that ∂U ∂λ varies by up to 20 kcal/mol and ∆∆G values by up to 10 kcal/mol across replicas for the chosen cases. Given the wide spread of these distributions, it is essential to perform ensemble simulations (often > 5 replicas) in order to control the aleatoric uncertainty associated with the ∆∆G estimations making the results reproducible. This can be demonstrated more clearly with the data in Table where results using a single replica have been compared with those using ensemble simulations for the entire dataset as well as various subsets of it. The accuracy of ensemble simulations is substantially better than any single replica, irrespective of the dataset taken.