Two-dimensional infrared (2D IR) spectroscopy is a modern experimental method used to sensitively probe molecular structures and dynamics. Since its inception by Hamm, Lim, and Hochstrasser 25 years ago, it has flourished into a broadly applicable experimental method. Applications in chemical biology include the elucidation of the structure and dynamics of proteins in solution. Recent experimental advances allow for studying biomolecules at interfaces, enable transient spectroscopy for following non-equilibrium protein dynamics, and implement spectro-electrochemistry to study redox-active proteins. Moreover, 2D IR spectroscopy is employed to investigate ground-state and excited-state dynamics in model systems such as metal carbonyl complexes and can be applied in material science to shine light on the ion dynamics in lithium ion batteries. However, the interpretation of 2D IR spectra critically depends on computations. Only by comparison to computational predictions, detailed information on the structure and dynamics can be extracted from experimental 2D IR spectra. Simply speaking, in 2D IR spectroscopy one performs a pump-probe experiment, in which a pump pulse excites a molecule from the vibrational ground state to a singly-excited vibrational state, followed by a probe pulse that results in stimulated emission or in an excitation to a doubly excited vibrational state. This way, 2D IR spectroscopy reveals the anharmonicities of the vibrational modes as well as couplings between them. By spreading out the spectroscopic information in two dimensions, 2D IR spectroscopy can reveal information that is hidden in conventional vibrational spectra. Furthermore, by varying the delay between the pump and the probe pulses, dynamical information can be assessed. Because it specifically probes anharmonicity, for a computational treatment of 2D IR spectroscopy one has to go beyond the harmonic approximation that is commonly used in theoretical vibrational spectroscopy. Common protocols for the simulation of 2D IR spectroscopy rely on empirical frequency maps in combination with vibrational exciton models (for an exhaustive review, see Ref. 34). Such models parametrize the frequencies and transition dipole moments of local vibrational modes as well as the harmonic couplings between them as functions of the molecular structure and the local electric field. Furthermore, the anharmonicity of the local modes is included as an empirical shift. This allows for the efficient calculation of the vibrational energy levels, which can be combined with molecular dynamics simulations and used as input for the simulation of 2D IR spectra. Such a protocol is efficient and has been successfully applied in numerous studies over the past 25 years. For commonly studied systems, reliable frequency maps that have been carefully parametrized with the help of quantum-chemical calculations are available, most importantly for the amide I vibrations in proteins and for the O-H and O-D stretching vibrations in water. However, the construction of frequency maps for new systems can be tedious and its parametrization requires a significant effort. Therefore, the reliance on such parametrizations has slowed down the progress in 2D IR spectroscopy.