Classical MD simulations were performed in three types of boxes: one identical to that used in AIMD simulations (i.e., L = 15 Å) and two larger boxes: L = 30 Å and L = 60 Å (Figure ). We found no significant differences in the solvation structure when comparing the 15 Å and 30 Å boxes (Figure ). Solvation free energy calculations of the charged methythiolate molecule were performed in both a 30 Å and 60 Å box. Importantly, we found no change in the calculated free energy between the two box sizes (Figure ) which suggests the 30 Å box is sufficiently large to minimize the non-neutral simulation cell artefact associated with the calculation. GROMACS 2023 was used to run all simulations. Simulations were carried out in the NVT ensemble with a constant temperature of 298 K, maintained using a Nose-Hoover thermostat with 3 ps coupling time. In the 30 Å and 60 Å boxes, long-range electrostatic interactions were calculated using the Particle-mesh Ewald method with a real-space cut-off of 1.2 nm and grid spacing of 0.12 nm with CHARMM and a real-space cut-off of 1.0 nm and grid spacing of 0.125 nm with Amber. For CHARMM the Lennard-Jones interactions were force-switched off between 1.0 and 1.2 nm, while for Amber, a cut-off at 1.0 nm was used and a dispersion correction was applied to the energy and pressure.