As summarized in Table , setting θ = 0.5 and p = 5, reduced the computational time for a TPA calculation by approximately one hour without compromising accuracy. Therefore, these settings were adopted for all subsequent calculations. Regarding the approximations on the QM/MM coupling, we examined the performance of single-center and ESPF methods. As shown in Table , ESPF with p = 0 and R exact = 10 provides the most effective balance between computational cost and accuracy on excitation energies, oscillator strengths, and TPA cross-sections. Based on these findings, the QM/MM coupling in this work was treated with ESPF with the aforementioned parameters. We calculated embedding potentials for both mfcc and cp3 approaches. Thereafter, excitation energies, oscillator strengths, and TPA cross-sections are obtained with PE-TD-DFT and CAM-B3LYP with Dalton program and the PE library. We included EEF effects and appropriate approximations (both for FMM and QM/MM couplings). The covalent bond crossing the quantum-classical border is treated as follows: The C α -C β bond is cut (Fig. ), and the dangling bond is saturated with a hydrogen atom (hereafter referred to as the hydrogen link atom). All charges in the environment that are within 0.8 Å of a nucleus in the quantum subsystem are distributed to the three nearest sites in the environment, while all other parameters are removed. To minimize the effect of ESO, ECPs are placed on atoms as described in Section 3.2. For the PE-TD-DFT calculation, several basis sets were studied, namely, 6-31+G*, cc-pVDZ, aug-cc-pVDZ, pcseg-1, aug-pcseg1, pcseg-2 and augpcseg-2, for both mfcc and cp3 approaches. The first excitation energy and two-photon absorption cross-section differences are shown, respectively, in Fig. and Fig. , where the reference is aug-pcseg-2. Timings for each basis set are shown in Fig. . We conclude that aug-pcseg-1 will be the basis set used thereafter since it is the computationally most effective and close to aug-pcseg-2 without compromising accuracy.