For the modeling of solvatochromism with an explicit representation of the solvent molecules, the quality of preceding molecular dynamics simulations is crucial. Therefore, the possibility to apply force fields which are derived with as little empiricism as possible seems desirable. Such an approach is tested here by exploiting the sensitive solvatochromism of p-nitroaniline, and the use of reliable excitation energies based on approximate second-order coupled cluster results within a polarizable embedding scheme. The quality of the various MD settings for four different solvents, water, methanol, ethanol, and dichloromethane, is assessed. In general, good agreement with the experiment is observed when polarizable force fields and special treatment of hydrogen bonding are applied.