where n is the total number of snapshots, σ 2 and µ are, respectively, the variance and mean of the TPA cross-sections σ TPA t at time t, and k is the lag (here k = 1). The TPA crosssections were calculated via PE-TD-DFT with CAM-B3LYP and 6-31+G * basis set, employing the cp3 embedding potential. This analysis was performed on the final 5 ps of the QM/MM MD trajectory, where 1000 snapshots were collected at 5 fs intervals. Fig. , shows the autocorrelation function fitted to an exponential decay model C(t) = A exp(-t/τ ), where A is the amplitude at t = 0 and τ is the characteristic decay time. The rapid decay of C(t) indicates that uncorrelated TPA structures can be found every 200 fs. To confirm this hypothesis and further evaluate the impact of snapshot selection on the spectra, we produced TPA spectra using 25, 50, 100 and 200 snapshots, corresponding to spacings of 1600, 800, 400 and 200 fs, respectively. While the positions and intensities are hardly affected by the number of snapshots (see Fig. ), the distribution of excitation energies, oscillator strengths and TPA cross-sections become more Gaussian-like with an increasing number of snapshots (Fig. ). For the generation of the MPA spectra, snapshots were extracted every 200 fs over the last 40ps QM/MM MD trajectory, resulting in a total of 200 snapshots. For each snapshot, we generate a cp3 embedding potential and conduct a PE-TD-DFT calculation with aug-pcseg-1 basis set. For each snapshot, this calculation is performed twice, where the core region is either the LYR-472 (QM-sampled) or LYR-225 (MM-sampled) moieties, as highlighted in Fig. .