Channelrhodopsin-2 (ChR2) is a light-gated ion channel widely used in optogenetics, a technique that enables precise control of neuronal activity by genetically engineering light-sensitive proteins into cell membranes. This protein exists in dimeric form, with each monomer containing a retinal Schiff base (RSB) moiety covalently bonded that undergoes trans-cis isomerization upon light absorption. However, the limited penetration depth of visible light in biological tissues motivates the use of multiphoton-absorption techniques, which enhance tissue penetration, improve focality, and reduce phototoxicity, thereby offering a promising alternative for optogenetic applications. In this paper, we present a fully atomistic multiscale methodology for computing the one-, two-, and three-photon absorption spectra of ChR2, where the protein, lipid bilayer, and solvent are explicitly considered throughout the workflow. This methodology integrates classical molecular mechanics (MM) molecular dynamics (MD), quantum mechanics/molecular mechanics (QM/MM)-MD, and fragment-based polarizable embedding (PE) to derive environment-specific PE potentials from the explicit protein-lipid-solvent environment. The final step in the methodology is to use these potentials to compute accurate spectra via PE-time-dependent density functional theory (PE-TD-DFT). Validation against experimental one-photon absorption spectra demonstrates excellent agreement. For the first time, we report the theoretical two- and three-photon absorption in ChR2, albeit without direct experimental comparison. We compare the multiphoton absorption (MPA) spectra where the two RSB moieties are sampled using classical MD and QM/MM-MD, respectively. The resulting spectral differences are attributed to variations in key structural parameters that we analyze and document.