## 2 Methodology
### Calculation of Energies and Spectral Density of FMO Complex  
In order to shed more light on the site energies and spectral density of the FMO complex, we recall the theoretical method used in [23], where a combination of quantum mechanical and molecular mechanics techniques have been used for computing them. Assuming that the quantities \(\mu_{k}\) and \(\mu_{m}\) are the strengths of the transition dipole moments for the \(k\)th and \(m\)th pigments, the coupling between the pigments can be expressed using the point dipole approximation as [23; 108; 109]:  
\[V_{km}=f\frac{\mu_{k}\mu_{m}}{R_{km}^{3}}[\vec{\mu}_{k}\cdot\vec{\mu}_{m}-3( \vec{\mu}_{k}\cdot\vec{\mu}_{km})(\vec{\mu}_{m}\cdot\vec{\mu}_{km})], \tag{12}\]where \(f\) is the protein screening factor and \(R_{ij}\) is the distance between transition dipoles with the direction \(\vec{n}_{km}\). For quantifying the spectral density, the autocorrelation function for each pigment can be computed from:  
\[C_{k}(t_{i})=\frac{1}{M}\sum_{l=1}^{M}\Big{[}\frac{1}{N-i}\sum_{j=1}^{N-i} \Delta E_{k,l}(t_{i}+t_{j})\Delta E_{k,l}(t_{j})\Big{]}, \tag{13}\]  
where the quantity \(\Delta E_{k,l}\) is the energy difference at time step \(t_{i}\) for the BChl \(k\) in monomer \(l\), which is given by:  
\[\Delta E_{k,l}=E_{k,l}(t_{i})-\big{\langle}E_{k,l}\big{\rangle}, \tag{14}\]  
where \(E_{k,l}(t_{i})\) is the instantaneous site energy and \(\big{\langle}E_{k,l}\big{\rangle}\) is the mean value of site energy. Furthermore, the autocorrelation function (\(C(t)\)) is related to the spectral density (\(J(\omega)\)) via the fluctuation-dissipation theorem as:  
\[C(t)=\frac{1}{\pi}\int_{\infty}^{\infty}d\omega\frac{e^{-i\omega t}J(\omega)} {1-e^{-\beta\omega}}, \tag{15}\]  
where \(J(\omega)\) is the spectral density and \(\beta\) is the inverse temperature denoted by \(\beta=1/(k_{B}T)\). Importantly, there are two different kinds of the models that can be used for spectral density (\(J(\omega)\)). One of them is the Drude model that can be expressed as [110]:  
\[J(\omega)=\frac{2\gamma\lambda\omega}{\omega^{2}+\gamma^{2}}, \tag{16}\]  
where \(\lambda\) denotes the reorganization energy that corresponds to the system-bath coupling strengths and \(\gamma\) is the cut-off frequency. The correlation function, approximated by a biexponential function with the Markovian white noise residue (WNR) ansatz, can be expressed as [110]:  
\[C(t>0)\approx ce^{-\gamma t}+c^{\prime}e^{-\gamma^{\prime}t}+2\Delta\delta(t). \tag{17}\]  
The associated parameters can be obtained, as in [110], by exploiting the Laurent expansion and [1/1] Pade approximation as:  
\[c=\frac{2\lambda}{\beta}\bigg{[}1-10(\beta\hbar\gamma/20)^{2}-\frac{2.45}{( \gamma^{\prime}/\gamma)^{2}-1}\bigg{]}-i\lambda\hbar\gamma, \tag{18}\]  
\[c^{\prime}=\frac{4.9\lambda\gamma\gamma^{\prime}}{\beta(\gamma^{\prime 2}- \gamma^{2})};\Delta=\beta\lambda\hbar\gamma/20. \tag{19}\]  
Equation (17) combined with Equations (18) and (19) is used for fitting the correlation function and spectral density can be computed by substituting \(\gamma,\gamma^{\prime}\) and \(\lambda\) into Equation (16). If the nuclear motion of the BChl system is considered in the protein environment, then the excitation energy transfer can be described by normal coordinates and the displacement of the nucleus can be described by [80]: \(\delta e^{\alpha}=\sum_{i}c_{i}^{\alpha}Q_{i}^{\alpha}\), where \(c_{i}^{\alpha}\) and \(Q_{i}^{\alpha}\) are the gradient of nuclear displacement and normal coordinates of the system, respectively. The spectral density of the system of the BChl can be computed by the determination of the correlation function. For providing better clarity, the computed results of BChl spectral density in [80] have been presented in Figures 3 and 4 utilizing various computational techniques. The spectra for the different kinds of BChls from the FMO complex for their mutated systems have been reported in [45].