Modeling the linear and nonlinear optical properties of chromophores in the condensed phase is highly challenging, as both the coupling of electronic excitations to nuclear vibrations and the interactions between the system and its complex environment have to be accounted for. A prominent example of this phenomenon are pigment-protein complexes, where the protein environment significantly alters the optical properties of embedded chromophores. Other examples include solvatochromic dyes, whose absorption and fluorescence spectra can undergo large changes depending on the solvent environment. In this talk, I will present a number of computational approaches designed to model these systems from first principles. The approaches are all based on large-scale electronic structure calculations including hundreds of atoms, to describe the chromophore and its polarizing environment on the same footing. It will be demonstrated that approximate treatments of the environment, such as through commonly used classical point charge models, can alter the couplings of electronic excitations to nuclear vibrations, with significant effects on predicted linear absorption spectra and nonlinear spectroscopy methods like 2D electronic spectroscopy (2DES). The first principles methods outlined in this talk are computationally expensive, generally requiring the calculations of thousands to tens of thousands of excitation energies along a molecular dynamics trajectory of the system in its complex environment. However, the computational cost can be significantly reduced by using machine learning (ML) approaches to predict optical excitation energies from test sets containing a few hundreds of structures, making the methods affordable enough for routine modeling of systems in complex condensed phase environments.