The excited state dynamics of [Fe(CN) 4 (bipy)] 2-is carried out using trajectory surface hopping simulations using SHARC and a vibronic coupling/molecular mechanics (VC/MM) model including all linear and selected quadratic coupling terms. The VC model has been parameterized 40 in a diabatic basis of 21 singlet and 20 triplet states with time-dependent density functional theory (TDDFT), where energies and gradients are obtained from GAUSSIAN16 and spin-orbit couplings from ORCA5. We employed the B3LYP* functional, using the def2-TZVP (for Fe) and def2-SVP all-electron basis sets (for other atoms), the GD3BJ empirical dispersion correction and the Tamm-Dancoff approximation. The charge transfer character is computed with TheoDORE, where-due to strong bonding-we consider the Fe(CN) 4 unit as the metal (M) and the bipy unit as the ligand (L). Further simulation details can be found in Sections S1.1 and S1.2 of the Electronic Supporting Information. Due to the high charge and polarity of [Fe(CN) 4 (bipy)] 2-, gasphase calculations have significant convergence problems. However, parametrizing the VC/MM model with implicit solvation would lead to double-counting of solvent shifts once the explicit solvent molecules are included. Hence, the VC parameters were generated using implicit solvation with a dielectric constant set to ε r = 1.77 (i.e., the square of the refractive index of water), such that the VC parameters take into account the polarizability of water solvent molecules due to their electron densities. The remaining polarizability of water is accounted for by the motion of the explicit water molecules within the VC/MM simulations (Section S1.1 and Figures ). Further, we include a selected subset of quadratic vibronic coupling parameters, i.e., state-specific frequency shifts in the form of γ (ii) kk parameters (Section S1.2 and Figure ) of normal modes that affect the equatorial Fe-ligand distances for all 3 MC states. Extending the VC model with these γ (ii) kk parameters significantly improved the accuracy of the MC states in preliminary gas-phase simulations. Both aspects, including implicit solvation into the parameters and including quadratic parameters, go beyond previous LVC/MM simulations. In order to yield enough trajectories to resolve the 3D solvation structure, we generated 30, 000 initial conditions of [Fe(CN) 4 (bipy)] 2-in water obtained from classical molecular dynamics simulations that used a flexible water model, as described elsewhere (Section S1.3). Out of the 30, 000 molecular dynamics snapshots, 4473 initial conditions were stochastically selected in a window between 487 nm to 497 nm (around 2.52 eV) centered on the peak of the first absorption band (Figure ). From them, 4366 initial conditions (97.6%) were excited into the adiabatic S 3 state, which has MLCT character. All 4473 initial conditions were propagated for 5000 fs on the lowest-lying six singlet and seven triplet adiabatic electronic states, which are enough for describing the excited state dynamics of the system. Details on the surface hopping settings can be found in Section S1.4.