from ∆G obtained from biased simulations using Equation S12. If it were possible to simulate a 'theoretically' infinitely dilute solution (the simulation box is big enough with adequate number of water molecules and La, Cl particles to reach bulk behavior), g ∞ (r) would converge to 1 and hence ϕ(r) would converge to 0 at r → ∞, which is the dissociated state in the ideal case. However, due to extreme long-range nature of Coulomb interactions and limited number of La, Cl particles, this is difficult to achieve even in simulation boxes with edge lengths of ∼30 Å . Hence, the dissociated state is taken to be the limiting continuum behavior where charged particles interact via screened Coulomb interaction given by ϕ analytical (r) = -q 2 /4πϵ 0 ϵ r r. FES obtained from simulation boxes of different sizes are rigidly shifted so that their tails coincide with the analytical solution at long distances, as shown in Figure . The difference in the log K ∞ eq values arise due to the difference in length of the tail coinciding with the analytical solution. However, once the FES is properly aligned to the analytical solution (indicating efficient sampling of all the metastable states), dependence of log K ∞ eq on cutoff distance is weak beyond ∼5 Å (corresponding to SShIP) at high temperatures, as shown in Figure . Therefore, results obtained from small simulation boxes used in AIMD simulations are acceptable at these conditions. The effect of box size may become important at lower temperatures where SSIPs are stabilized in addition to CIPs and SShIPs. Since SSIPs are formed at ∼7.5 Å-8 Å efficient sampling of these ion pairs is no longer possible in small simulation boxes.