QM model calculations. All QM model calculations were performed with the Gaussian 16 software. The geometries of interested species were fully optimized in conjunction with the SMD continuum solvation model at the B3LYP-D3/def2-SVP level of theory. The energies were further refined with the larger basis set def2-TZVPP for all atoms. Harmonic frequency calculations were performed at the same level of theory as the optimizations in order to estimate the zero-point energies as well as the thermal and entropic corrections. In addition, we also tested the 𝜔B97XD functional for calculation of the electronic coupling. With 𝜔B97XD, the calculated electronic coupling for the second ET step is 0.015 eV, which is close to value of 0.021 eV obtained by B3LYP. As QM-optimization may lead to shortening of the donor-acceptor distance compared to the one within protein environment, we have adjusted the electronic couplings based on the difference between the QMoptimized ET distance and the averaged ET distance from the QM/MM-MD simulation (Eq.2 and Eq. 5 in SI).