  The $GW$-Bethe-Salpeter Equation (BSE) method is promising for calculating
the low-lying excited states of molecular systems. So far, it has only been
applied to rather small molecules, and in the commonly implemented diagonal
approximations to the electronic self-energy it depends on a mean-field
starting point. We describe here an implementation of the self-consistent and
starting-point independent quasiparticle self-consistent (qs$GW$)-BSE approach
which is suitable for calculations on large molecules. We herein show that
eigenvalue-only self-consistency leads to an unfaithful description of certain
excitonic states for Chlorophyll dimers while the qs$GW$-BSE vertical
excitation energies (VEE) are in excellent agreement with spectroscopic
experiments for Chlorophyll monomers and dimers measured in the gas phase. On
the other hand, VEEs from time-dependent density functional theory calculations
tend to disagree with experimental values and using different range-separated
hybrid (RSH) kernels changes the VEEs by up to 0.5 eV. We use the new
qs$GW$-BSE implementation to calculate the lowest excitation energies of the
six chromophores of the photosystem II (PSII) reaction center (RC) with nearly
2000 correlated electrons. Using more than 11000 (6000) basis functions, the
calculation could be completed in less than 5 (2) days one a single modern
compute node. In agreement with previous TD-DFT calculations using RSH kernels
on models that do also not include environment effects, our qs$GW$-BSE
calculations only yield states with local character in the low-energy spectrum
of the hexameric complex. Earlier work with RSH kernels has demonstrated that
the protein environment facilitates the experimentally observed
interchromophoric charge transfer. Therefore, future research will need to
combine correlation effects beyond TD-DFT with an explicit treatment of
environment electrostatics.
