is more stable by 0.580 eV. However, the picture changes dramatically when one adds higher-body terms to the expression of the total interaction energy: The HS state becomes the more stable state by 0.161 eV. Previous solvation studies have used MM and electrostatic QM/MM models. The former of these models neglects many-body effects such as induced polarization everywhere, and the latter only in the MM region. these models produce solvation structures that are generally less strongly bound than what is predicted by fully QM AIMD simulations, which includes many-body interaction. Based on the analysis above, we propose that some of this discrepancy is caused by many-body effects within the solvation shell of the complex. Interestingly, as the QM/MM MD simulations did also not fully reproduce the AIMD results, further testing this hypothesis using polarizable embedding QM/MM models could prove insightful. For this specific complex, the aforementioned role of the pose-A water molecules in their close vicinity to the iron center could mean that simply including these three molecules in the QM subsystem might already improve the model significantly.