We present a general and comprehensive theoretical and computational framework for modeling ultrafast multidimensional infrared spectra of a vibrational excitonic system in liquid solution. Within this framework, we describe the dynamics of the system in terms of a quantum master equation that can account for population relaxation, dephasing, coherence-to-coherence transfer, and coherence-to-population transfer. A unique feature of our approach is that, in principle, it does not rely on any adjustable fitting parameters. More specifically, the anharmonic vibrational Hamiltonian is derived from ab initio electronic structure theory, and the system bath coupling is expressed explicitly in terms of liquid degrees of freedom whose dynamics can be obtained via molecular dynamics simulations. The applicability of the new approach is demonstrated by employing it to model the recently observed signatures of coherence transfer in the two-dimensional spectra of dimanganese decacarbonyl in liquid cyclohexane. The results agree well with experiment and shed new light on the nature of the molecular interactions and dynamics underlying the spectra and the interplay between dark and bright states, their level of degeneracy, and the nature of their interactions with the solvent.