The soft modes, including the flexible degrees of freedom (DoFs) of both the solute and solvent, are sampled using classical MD. In contrast, the stiff modes are treated at the vibronic level using harmonic potentials specific to each configuration of the soft modes. To this end, analytic formulations within the BO approximation of the vibronic spectrum have been derived, and are implemented in our open source code F C classes3. The final Ad-MD|gVH spectrum is derived by averaging signals computed at the QM level for the stiff coordinates across the configurations obtained from the MD trajectory, which explicitly incorporates the equilibrated solvent. Further details on the methodology and related computations are given in Section S4