The current consensus in the field of molecular dynamics simulation is that increasing the length of a single simulation leads to improvement in the accuracy and precision of calculated expectation values. On the basis of chaos theory, and the fact that the ergodic theorem cannot hold for molecular dynamics simulations on accessible timescales, we probed this assumption and provided direct evidence that individual trajectories do not suffice for deriving precise, reproducible and accurate results for protein-ligand systems. We showed on the contrary that ensembles are essential for the calculation of statistically robust results, regardless of the length of simulation for this class of systems. On comparing the proteinligand contact frequency distributions from ten or twelve independent 10 µs trajectories, 90% or more pairs of trajectories had significantly different distributions of ligand-protein residue interactions. The principles and findings of this study are not just confined to ligand-protein systems and free energy calculations, being based on the chaotic nature of any dynamical system which displays an equilibrium state are more widely applicable to molecular dynamics in general and hence should be accounted for in all MD based applications regardless of the particular domain of interest. To investigate the effect of this uncertainty on the value of a one dimensional macroscopic observable, we analysed the same set of trajectories in order to determine ligand binding free energies and their associated statistical distributions. The specific method which we used for ligand binding free energy calculations was taken from Pan et al. In their paper, the authors reported strong correlation to FEP calculations but poor correlation to experiment, stating this poor correlation may be attributable to force field inaccuracies. In the present study, we demonstrated that separate trajectories lead to the computation of completely different results, differing by up to 7.26 kcal/mol. Our study conclusively demonstrates that binding free energies from individual simulations are inherently non-precise, non-reproducible and do not yield chemical accuracy (± 1 kcal/mol ). Clearly, long-timescale trajectories probe an insufficient number of microstates to effectively sample the phase space. In turn, the lack of agreement with experiment should not necessarily be attributed to force-field inaccuracies. This is a paramount example of the importance of taking aleatoric uncertainty fully into account.