The Fe complex was parametrized using MCPB.py, as described elsewhere, although its intra-molecular potential terms were not used (vide infra). Partial charges for each electronic state were generated using the CHELPG model, shown previously to be one of the bestperforming methods for these systems. The complex was then solvated in a box of 7329 water molecules modelled with the TIP4PEW potential, using the 'leap' program from the 'Ambertools' toolbox, which also added two Cl -counterions, to charge-neutralize the system. To be consistent with previous methodology, focusing on solvation-shell geometries rather than internal solute geometries/dynamics, the complex was restrained to its relaxed B3LYP-D3(BJ)/def2-TZVP geometry (for each electronic state) using 500 kJ/mol positional restraints on each atom in the complex, and the counterions were restrained 21-23 Å away from the center of the complex to avoid them affecting the inner solvation shell.