In general, spontaneous ligand binding methods work by initiating a molecular dynamics protein-ligand system from a configuration where the ligand is placed at some distance from the surface of the protein. During the simulation, the ligand explores the surface of the protein and binds to potentially druggable sites which may be orthosteric, allosteric or even cryptic in nature. By analysing the trajectories using methods such as Markov state model (MSM) analysis, thermodynamic and kinetic observables which are of key importance to the process of drug discovery can be extracted from the data. These include binding free energies, dissociation constants (K d ) and on and off rates of binding (K on and K of f ). When conducting these studies, the question arises as to whether the trajectory has sufficiently sampled phase space such that the probability distribution of the trajectory has converged to the equilibrium probability distribution of the protein-ligand system. Only this distribution would allow the true expectation values of the observables to be obtained. To sample the phase space, one of two distinct approaches is usually followed. In the first, which we term the "long-timescale" regime, authors report several microsecond timescale simulations and from these, compute the observables of interest. These observables include 3-dimensional ligand occupancy maps, ligand binding free energies, along with ligand association mechanisms and pathways. In the second regime, termed "adaptivesampling", many simulations of shorter timescale are executed, and new simulations are adaptively initiated from specific simulation snapshots in order to "more thoroughly" explore regions of phase space that are of interest. Incidentally, many studies from the second regime report aggregate simulation times that fall in the microsecond timescale; this is misleading as performing a single simulation of that duration is not the same as we will discuss in detail in the current study. We note that some studies also combine the two techniques, using adaptive sampling to initiate new, "short" simulations from long-timescale simulations that are stuck within non-productive kinetic traps. Generally, this approach is taken in order to converge transition probabilities between metastable states that are identified during Markov state modelling. We would like to point out here that there are several accelerated sampling protocols that involve performing "ensembles" of simulations. These include methods that do not employ any external force or heating, and just enhance sampling by performing multiple independent MD simulations concurrently with different starting conditions. Examples include ensemble dynamics, Markov state model (MSM), weighted ensemble (WE) and multilevel splitting (MS) methods. Although these methods involve performing "replicas" and generating "ensembles", the fundamental question is whether we get the same answer (within error bars) on repeating the entire protocol using one of the above methods. Given that the dynamics is chaotic, it is expected that this is not the case and ensembles must be used as each execution of such a protocol would have a different initial condition. One example are replica exchange methods that also involve performing multiple MD simulations in parallel (so called "replicas"). We have shown in previous work that on repeating a replica exchange calculation multiple times, we indeed observe variation in the outcome and hence it is necessary to perform ensembles of the entire protocol (which itself contains "replicas") to perform a systematic uncertainty quantification (UQ). Similar studies are required for other methods involving "ensembles" in order to properly assess UQ in those cases.