In detail, the bound vanish stage time series generally showed no substantial drift over the final 20 ns (from 10 -30 ns, Figures and). Therefore, the "true" infinite sampling time free energy changes were calculated as the mean over the last 20 ns, and were subtracted from the time series. The initial transients were generally well-fit by a single exponential decay (Figure ), although some quickly-decaying bias was evident for most systems at short timescales. To account for this, a second exponential with short half-life was fit after subtracting the first exponential, and retained only if the pre-exponential term was positive 13 (Figure ). Autocovariance functions were calculated for each system using the final 20 ns portions of the time series. The first two points were calculated directly from the averaged time series, while the remaining points were obtained by interpolating the convex gamma function computed following Geyer (Figures and). Figure shows that the final 20 ns portions of the time series were well-modelled by Gaussian distributions (although some systems showed small but significant deviations from normality). Therefore, synthetic ensembles of trajectories were generated for each system starting from vectors of uncorrelated standard normal noise. 1000 noise vectors were generated for each system. The Cholesky decompositions of the autocovariance matrices were obtained and the resulting lower-triangular matrices were used to reintroduce the desired correlation structure to the uncorrelated noise vectors. The exponentials were then added to introduce the initial transients. Examples of synthetic time series are shown in Figures and.