MD simulations are performed in the N V T ensemble, that is, with a constant number of particles, N , volume, V and temperature, T , and with a time step of 0.5 fs. After initial equilibration, the simulations are run for approximately 5 ns. To compare the fitted PIM potential against available non-polarizable force fields, MD simulations are carried out with the same boxes using the pair potential reported in Migliorati et al. and Smith and Dang , referred to as N-POL in the following. Details of the N-POL potential are described in Table . These potentials have been shown to accurately predict the structural and thermodynamic properties of aqueous La 3+ at ambient conditions . To check the transferability of the PIM potential to subcritical conditions, additional MD simulations are performed at 473 K and 400 bar. This temperature and pressure are chosen to ensure solvent densities similar to the training set and to facilitate direct comparison with AIMD data from Guan et al. .