Here G is the Gibbs free energy and ∆G is the change in Gibbs free energy between two states A and B. ∆G is calculated by integrating ∂u(λ,x) ∂λ λ over the range of λ and this integration is performed numerically. It is worth highlighting that equation 2 is only strictly valid in the thermodynamic limit, when both left-hand-side and right-hand-side terms are unique numbers with no fluctuations. However, practically speaking, we work with finite systems and sample only a fraction of the full conformational space, which makes these quantities stochastic variables . This implies that both the free energy as well as its derivative will have a distribution of values. Therefore, it is necessary to get the expectation value of these quantities using ensemble methods. The brackets . λ denote an ensemble average in a thermodynamic state defined by the value of λ. To compute this ensemble average the configurations of particles can be sampled using Monte Carlo or MD methods and from these sampled configurations values of the potential are calculated and averaged. The traditional approach has been to perform a single "long" MD simulation to proxy ensemble averaging with time averaging. However, as discussed already, this is not reliable due to the extreme sensitivity of results obtained, arising from the initial conditions which are controlled by the random seeds used to initiate simulations. Thus, ensemble simulations are required to generate the ensemble of conformations in order to estimate an average and distribution of the calculated ∆G. In this work, the same idea of performing ensemble simulation to get the expectation value of the distribution of ∆G by performing stochastic integration of the distributions of ∂u(λ,x) ∂λ λ is applied using the TIES methodology.