Using classical molecular dynamics simulations, the 2D infrared (IR) spectroscopy of CN⁻ solvated in D₂O is investigated. Depending on the force field parametrizations, most of which are based on multipolar interactions for the CN⁻ molecule, the frequency-frequency correlation function and observables computed from it differ. Most notably, models based on multipoles for CN⁻ and TIP3P for water yield quantitatively correct results when compared with experiments. Furthermore, the recent finding that T ₁ times are sensitive to the van der Waals ranges on the CN⁻ is confirmed in the present study. For the linear IR spectrum, the best model reproduces the full widths at half maximum almost quantitatively (13.0 cm⁻¹ vs. 14.9 cm⁻¹) if the rotational contribution to the linewidth is included. Without the rotational contribution, the lines are too narrow by about a factor of two, which agrees with Raman and IR experiments. The computed and experimental tilt angles (or nodal slopes) α as a function of the 2D IR waiting time compare favorably with the measured ones and the frequency fluctuation correlation function is invariably found to contain three time scales: a sub-ps, 1 ps, and one on the 10-ps time scale. These time scales are discussed in terms of the structural dynamics of the surrounding solvent and it is found that the longest time scale (≈10 ps) most likely corresponds to solvent exchange between the first and second solvation shell, in agreement with interpretations from nuclear magnetic resonance measurements.