The umbrella sampling method was used to compute the potential of mean force for the aggregation between pairs of NPs in both ILs. On this method, a bias is introduced to force the sampling of the desired reaction coordinate, in this case, the distance r between the centers of mass of the two NPs. To sample all the distances between the contact distance and the maximum one enabled by the size of the simulation box, a total of 215 simulations were performed for each pmf by changing the minimum r min of the harmonic potential used as the bias potential (Equation ) from 4.40 to 12.35 nm. Hence, those calculations reproduced the dissociation between the NPs starting from the structure with the lowest energy found for the dimer. Each simulation with the slightest high r min started from the final structure of the previous simulation, so the change in the potential corresponds only to a very small perturbation on the system, and was performed for 44 ns for r min < 6.0 and for 50 ns for r min ≥ 6.0, with the first 4 ns of each simulation discarded from analyses. The probability distributions obtained in each simulation were combined by the weighted histogram analysis method (WHAM) as implemented in Gromacs to generate the potential of mean force.